This analysis notebook accompanies the paper Momentary mood and arousal detect real-life prolonged stress by Tutunji et al. 2021, and utilizes data from the study “Stress Reslience and the Brain in Medical Students (STRAIN-MD)”. Study design details can be founed in the original article, but a brief summary is provided here to assist in reading comprehension. In this study, participants completed two weeks of Ecological Momentary Asessments (EMA/ESM) coupled with Ecological Physiological Assessments (EPA). One of these weeks culminated in a high-stakes exam, while the other was outside of examination periods thus providing a stress and control week respectively. EMA consisted of questions regarding subjective stress, and positive and negative affect. EPA within the context of this study consisted of a wearable device that collected continuous measures of physiological arousal. Using these measures, we first validate our stress exposure by examining changes in subjective stress ratings between the two week. We then demonstrate the impact of stress exposure on mood and arousal measures. We then explore our findings by looking at continuous associations between subjective stress and outcomes. Finally, we demonstrate the feasibility of Mood and EPA measures in prediction of stress and control weeks using machine learning models.
All custom functions and code can be found in associated GitHub directory. Libraries used listed at the end of the document. In addition to the libraries, we also used prespecified themes for ggplots that make some plots easier to draw. These are all run in the set-up chunks.
# Libraries to import and clean data
library(readxl)
There were 50 or more warnings (use warnings() to see the first 50)
library (tidyverse)
library(plyr)
library(broom)
library (hablar)
library(janitor)
library(data.table)
library(zoo)
library(scales) # Rescaling made easy
# Library for EMA
library (remotes)
#remotes::install_github("wviechtb/esmpack") # commented, needs to be installed from github
library(esmpack)
library(boot)
# Basic Stats
library (psych)
library (psycho)
library(Hmisc)
library(RcmdrMisc)
library(corrr)
library(ppcor)
# Mixed Model
library (lmerTest)
library(optimx) # Optimizer for lm
library(olsrr) # Diagnostic for lmer
library(jtools) # Uses summ function for readable models that look pretty
library(sjPlot) # more modeling display
library(effectsize) # Get effect sizes from tstat
library(r2glmm) # Compute effect size directly
library(performance)
# Parallel processing
library(future)
library(doSNOW)
library(randomForest)
#Plotting
library(ggplot2) # Standard plotting
library(qgraph)
library(ggpubr)
library(knitr)
#PCA
library(factoextra)
#Fonts and Pretty things
library(stargazer)
library(extrafont)
#font_import()
loadfonts(device="pdf")
# Theme 1
ggtheme <- theme (text=element_text(size=16, family="Cambria"),
legend.position = "none",
axis.text.y = element_blank(),axis.text.x = element_text(size=16),
axis.ticks.y = element_blank(),
plot.title = element_text(size=16, hjust=0.5),
panel.background = element_rect(fill="transparent"),
panel.grid.minor.y = element_line(size=3),
panel.grid.major = element_line(colour = "aliceblue") )
# Theme 2
ptheme <- theme(text=element_text(size=11, family="Calibri"),
axis.line = element_line(size = 1, colour = "grey"),
panel.background = element_rect(fill="transparent"),
plot.background = element_rect(fill = "transparent", color = NA), # bg of the plot
panel.grid.minor.y = element_line(colour="grey95"),
panel.grid.major.y = element_line(colour = "grey95"))
# Set number of cores
NoCores <- as.numeric(availableCores())-1
# Load custom functions
source("functions.R")
We first load in the descriptive file here, and perform some cleaning to make the file more human friendly. While no identifying information is directly stored in the descriptive file, there still is a potential that associations can be made. This section is redacted from the public documentation for this reason.
# Redacted
Using python, we derived physiological arousal features in 10-minute time windows prior to each EMA survey from the raw EPA wristwatch data for each participant. Preprocessing is done through python with the pyphysio package which has some tools for processing bio-signals (refer to link for example of the standard processing pipeline). Several measures are derived at this point from the physiology files that scope electrodermal activity, heart rate, temperature, and movement. The data frame with the derived EPA measures, and the raw unporcessed EMA measures contain the following features:
Due to privacy issues, only a specific file (and not the full file) with anonymized categorical variables for the program are presented in this code. All data presented can be used for replication of the findings. Generated variables include different suffixes, which each represent an alternative mode of calculating the measures of interest. Below is a description of what each suffix stands for:
Below is the anonymized data frame used for the analysis:
#EMA_Data <- read.csv("data/EMA_Clean_Anon.csv")
print(EMA_Pub)
Before running any stats, we should take a look at some descriptive stats to get an overall feel of our data. We do this in two parts, looking at the population descriptives, and then the physiology descriptives.
Lets first check compliance rates to see how well our subjects adhered to the sampling scheme. We first check the compliance rates uncorrected for the time at which they were acquired, followed by correction of the time limits. We set this to 1 hour since we only have 6 surveys, and our subjects had varying time schedules.
Compliance rates uncorrected for time off-sets
# Recod compliance variable
EMA_Data$completed[EMA_Data$survey_progress==100] <- 1
Unknown or uninitialised column: `completed`.
# Per subject
calc.nomiss(completed, castor_record_id, data=EMA_Data, prop=T)
sub_001 sub_002 sub_003 sub_004 sub_005 sub_006 sub_007 sub_008 sub_009 sub_010 sub_011 sub_012 sub_013 sub_014 sub_015 sub_016
0.9759036 0.8902439 0.9259259 0.8148148 0.9259259 1.0000000 0.9518072 0.9135802 0.9880952 0.6125000 1.0000000 0.9756098 0.9880952 0.9756098 0.6835443 1.0000000
sub_017 sub_018 sub_019 sub_020 sub_021 sub_022 sub_023 sub_024 sub_025 sub_026 sub_027 sub_028 sub_029 sub_030 sub_031 sub_032
0.8421053 0.7307692 0.9759036 0.9113924 0.8734177 0.9113924 0.9759036 0.9879518 0.9487179 0.9878049 0.9500000 0.9753086 0.9240506 0.9480519 0.8271605 0.7887324
sub_033 sub_034 sub_035 sub_036 sub_037 sub_038 sub_039 sub_040 sub_043 sub_044 sub_045 sub_046 sub_047 sub_048 sub_050 sub_052
0.9210526 0.9240506 0.9382716 0.9375000 0.9743590 0.6562500 0.9878049 0.9506173 0.9638554 0.9518072 0.9166667 0.9634146 0.9382716 0.8795181 0.9350649 0.9500000
sub_053 sub_054 sub_057 sub_058 sub_059 sub_060 sub_061 sub_062 sub_063 sub_064 sub_065 sub_066 sub_067 sub_068 sub_069 sub_070
0.8974359 0.9390244 0.9375000 0.9506173 0.9629630 0.9500000 1.0000000 1.0000000 0.8947368 0.9634146 0.9880952 0.9759036 1.0000000 0.9761905 0.8028169 0.9367089
sub_071 sub_072 sub_073 sub_074 sub_075 sub_076 sub_077 sub_078 sub_081 sub_082 sub_083 sub_084 sub_085 sub_086 sub_087 sub_088
0.9625000 0.9125000 0.9350649 NA 0.9761905 0.9879518 0.9518072 1.0000000 0.9880952 0.8780488 0.8170732 0.8571429 NA NA NA NA
sub_089 sub_105 sub_115 sub_141 sub_142 sub_149 sub_151 sub_153 sub_155 sub_156 sub_161
NA 0.7230769 1.0000000 0.9880952 0.9629630 0.9200000 1.0000000 NA 0.8400000 0.9200000 NA
# more summary statistics for compliance
summary(calc.nomiss(completed, castor_record_id, data=EMA_Data, prop=TRUE))
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.6125 0.9119 0.9500 0.9249 0.9759 1.0000 8
Without correcting for time offsets, our EMA measures seem to be nice. Lets see if I correct for timing if it makes a difference next what will happen. We will give participant one hour to finish the survey. If they don’t finish it in that time window, we will count it as missing.
Compliance rates correcting for time of completion. First drop any survey that is not complete more than 80%. I then create a variable that tells as when the survey was completed in terms of time (not date). I then conditionally set the completion variable to zero if the survey does not fall within the completed time window given the beep.
# Set up some stuff for replace and refill
EMA_Data$survey_completed_on <- as.character(EMA_Data$survey_completed_on)
EMA_Data$survey_completed_time <- as.ITime(EMA_Data$survey_completed_on)
EMA_Data$completed[EMA_Data$survey_progress < 80] <- 0
# Replace and refill conditionals for completion, factoring times
EMA_Data$completed[EMA_Data$survey_completed_time >= (as.ITime("10:30:00")) & EMA_Data$ema_beep==1] <- 0
EMA_Data$completed[EMA_Data$survey_completed_time >= (as.ITime("13:00:00")) & EMA_Data$ema_beep==2] <- 0
EMA_Data$completed[EMA_Data$survey_completed_time >= (as.ITime("15:30:00")) & EMA_Data$ema_beep==3] <- 0
EMA_Data$completed[EMA_Data$survey_completed_time >= (as.ITime("18:00:00")) & EMA_Data$ema_beep==4] <- 0
EMA_Data$completed[EMA_Data$survey_completed_time >= (as.ITime("20:30:00")) & EMA_Data$ema_beep==5] <- 0
EMA_Data$completed[EMA_Data$survey_completed_time >= (as.ITime("23:00:00")) & EMA_Data$ema_beep==6] <- 0
# Reformat 0 to Nan, and get new completion rates
EMA_Data$completed[EMA_Data$completed == 0] <- NA
calc.nomiss(completed, castor_record_id, data=EMA_Data,prop=TRUE)
sub_001 sub_002 sub_003 sub_004 sub_005 sub_006 sub_007 sub_008 sub_009 sub_010 sub_011 sub_012 sub_013 sub_014 sub_015 sub_016
0.9638554 0.7073171 0.7160494 0.6913580 0.9259259 0.9761905 0.9277108 0.8395062 0.9404762 0.5625000 0.8641975 0.9512195 0.9404762 0.7682927 0.4050633 0.9404762
sub_017 sub_018 sub_019 sub_020 sub_021 sub_022 sub_023 sub_024 sub_025 sub_026 sub_027 sub_028 sub_029 sub_030 sub_031 sub_032
0.8421053 0.7051282 0.9277108 0.8734177 0.7721519 0.8101266 0.9638554 0.9879518 0.9487179 0.9756098 0.9375000 0.9506173 0.9240506 0.9480519 0.8271605 0.7887324
sub_033 sub_034 sub_035 sub_036 sub_037 sub_038 sub_039 sub_040 sub_043 sub_044 sub_045 sub_046 sub_047 sub_048 sub_050 sub_052
0.8815789 0.8354430 0.9259259 0.9250000 0.9615385 0.4687500 0.8658537 0.9382716 0.9518072 0.9397590 0.7976190 0.9390244 0.9382716 0.8795181 0.8961039 0.9500000
sub_053 sub_054 sub_057 sub_058 sub_059 sub_060 sub_061 sub_062 sub_063 sub_064 sub_065 sub_066 sub_067 sub_068 sub_069 sub_070
0.8717949 0.9268293 0.9375000 0.9506173 0.9259259 0.9375000 0.9523810 1.0000000 0.8815789 0.9146341 0.9166667 0.9397590 0.9285714 0.6428571 0.5492958 0.7721519
sub_071 sub_072 sub_073 sub_074 sub_075 sub_076 sub_077 sub_078 sub_081 sub_082 sub_083 sub_084 sub_085 sub_086 sub_087 sub_088
0.9125000 0.8500000 0.7922078 NA 0.9404762 0.8915663 0.9156627 0.8809524 0.9404762 0.8048780 0.7317073 0.8571429 NA NA NA NA
sub_089 sub_105 sub_115 sub_141 sub_142 sub_149 sub_151 sub_153 sub_155 sub_156 sub_161
NA 0.7076923 0.8780488 0.8452381 0.9629630 0.8400000 0.9642857 NA 0.6133333 0.9200000 NA
summary(calc.nomiss(completed, castor_record_id, data=EMA_Data,prop=TRUE))
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.4051 0.8313 0.9157 0.8638 0.9405 1.0000 8
We also check the completion rates for periods where we have detected usable physiology in our pipeline. We do this by estimating the number of times where we successfully detected skin conductance,
summary(calc.nomiss(sc_tonic_mean, castor_record_id, data=EMA_Data,prop=TRUE))
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.0000 0.7037 0.7711 0.7337 0.8373 0.9762 8
When we correct for the timing, we actually still have decent completion rates. Finally we filter out all incomplete surveys and those outside our time windows for our overall analysis.
EMA_Data <- EMA_Data %>% dplyr::filter(!is.na(completed))
We next take a look at population descriptive that we would like to report.
describe(EMA_descr$Sex)
EMA_descr$Sex
n missing distinct
83 0 2
Value Female Male
Frequency 51 32
Proportion 0.614 0.386
This is interesting for other analyses we want to do with other data
describe(factor(EMA_descr$Contraceptive_use))
factor(EMA_descr$Contraceptive_use)
n missing distinct
83 0 4
Value c Male No Yes
Frequency 8 32 19 24
Proportion 0.096 0.386 0.229 0.289
describe((EMA_descr$Program))
(EMA_descr$Program)
n missing distinct
83 0 2
Value BSc_BioMed BSc_Med
Frequency 22 61
Proportion 0.265 0.735
describe(EMA_descr$First_Week)
EMA_descr$First_Week
n missing distinct
83 0 2
Value Control-First Exam-First
Frequency 56 27
Proportion 0.675 0.325
# Number of days between weeks
summary(abs(EMA_descr$Column1))
Min. 1st Qu. Median Mean 3rd Qu. Max.
10.00 14.00 15.00 15.96 16.00 33.00
Asking whether the surveys disrupted the participants.
psych::describe(EMA_Data$physical_disruption)
ggplot(EMA_Data, aes(x=physical_disruption)) + geom_histogram(bins = 7)
After getting the population descriptive, we should next derive some E4 Quality measures. A lot of the data which is poor quality is already removed early on in the preprocessing pipeline. Instead we can check how many recordings were captured for the IBI data, and the overall movement during testing.
This is a measure of what duration of time the detected heart beats cover. So within a 10minute time window, we would want to detect 10 minutes of total IBI time if we had 100% signal. We see that for the mean is 27% of time with beats detected instead. This indicates that we do not have enough data for temporal domain analysis.
ggplot(EMA_Data, aes(x=ibi_based_quality, colour=ibi_based_quality)) + geom_histogram()
count(EMA_Data$ibi_based_quality > 0.3)
psych::describe(EMA_Data$ibi_based_quality)
Next we will check if the IBI quality is related to the ACC data. This gives us a good idea of how much motion is affecting our data. We see it has a big impact. When people are not moving, we are able to detect the entire IBI window. This means we should include movement in our models.
ggplot(EMA_Data, aes(x=ibi_based_quality, y=acc_delta)) + geom_smooth(method="lm", fullrange=F ) + coord_cartesian(ylim = c(0,10), xlim=c(0, 2.5)) + ggtheme + theme(axis.ticks.y=element_line(size=1), axis.text.y = element_text(size=12)) + ylab("ACC Delta") + xlab ("IBI Based Quality")
Next, we check if the IBI Quality effects the SD of the HR measures. It looks like it does have an impact, meaning we likely cant use that measure.
ggplot(EMA_Data, aes(x=ibi_based_quality, y=hr_sd)) + geom_smooth(method="lm", fullrange=F ) + ggtheme + theme(axis.ticks.y=element_line(size=1), axis.text.y = element_text(size=12)) + ylab("HR SD") + xlab ("IBI Based Quality")
summary(lm(ibi_based_quality ~ hr_sd, data=EMA_Data))
Call:
lm(formula = ibi_based_quality ~ hr_sd, data = EMA_Data)
Residuals:
Min 1Q Median 3Q Max
-0.42413 -0.16391 -0.05473 0.11399 1.45944
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.4722122 0.0061893 76.3 <2e-16 ***
hr_sd -0.0239518 0.0006238 -38.4 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.2216 on 4417 degrees of freedom
(1316 observations deleted due to missingness)
Multiple R-squared: 0.2502, Adjusted R-squared: 0.2501
F-statistic: 1474 on 1 and 4417 DF, p-value: < 2.2e-16
Next we will check if the skin conductance signals are related to the ACC. Basically this is just doing exploratory analysis to double check whether movement also impacts skin conductance measures. We see here that it also does indeed impact our measures.
fig.sc_tonic <- ggplot(EMA_Data, aes(x=sc_tonic_mean, y=acc_delta)) + geom_smooth(method="lm", fullrange=F ) + ggtheme + theme(axis.ticks.y=element_line(size=1), axis.text.y = element_text(size=12)) + ylab("ACC Delta") + xlab ("SC Tonic")
fig.sc_mag <- ggplot(EMA_Data, aes(x=sc_phasic_mag, y=acc_delta)) + geom_smooth(method="lm", fullrange=F ) + ggtheme + theme(axis.ticks.y=element_line(size=1)) + xlab ("SC Mag")
fig.sc_num <- ggplot(EMA_Data, aes(x=sc_phasic_num, y=acc_delta)) + geom_smooth(method="lm", fullrange=F ) + ggtheme + theme(axis.ticks.y=element_line(size=1)) + xlab ("SC Num")
fig.sc_auc <- ggplot(EMA_Data, aes(x=sc_phasic_auc, y=acc_delta)) + geom_smooth(method="lm", fullrange=F ) + ggtheme + theme(axis.ticks.y=element_line(size=1)) + xlab ("SC AUC")
fig.sc_dur <- ggplot(EMA_Data, aes(x=sc_phasic_dur, y=acc_delta)) + geom_smooth(method="lm", fullrange=F ) + ggtheme + theme(axis.ticks.y=element_line(size=1)) + xlab ("SC Dur")
fig.temp_mean <- ggplot(EMA_Data, aes(x=physical_excercise_dur, y=acc_delta, colour="red")) + geom_smooth(method="lm") + ggtheme
ggarrange(fig.sc_tonic, fig.sc_mag, fig.sc_num, fig.sc_auc, fig.sc_dur, fig.temp_mean)
fig.sc_mag <- ggplot(EMA_Data, aes(x=sc_phasic_mag, y=sc_phasic_num)) + geom_smooth(method="lm", fullrange=F ) + ggtheme + theme(axis.ticks.y=element_line(size=1)) + xlab ("SC Mag")
fig.sc_mag
We finally plot the graphs of the EMA items over time in each of the weeks for some exploratory analysis.We see here that there is an effect of time (higher stress in mid-day), and an effect of day. The closer they are to exams, the more stressed they are. This means we need to include these in our models.
vars.ema <- c( "event_tot_z", "activity_tot_z", "social_tot_z", "physical_tot_z", "mood_positive_z", "mood_negative_z")
for ( i in vars.ema) {
p1 <- ggplot(EMA_Data, aes(y=get(i), x=ema_survey, colour=Week_Type, fill=Week_Type)) + geom_line(stat="summary") + geom_ribbon(stat="summary", alpha=0.3, color=NA) + scale_x_continuous() + ylab(i)
print(p1)
}
.
After plotting our variables over all the week, it will be interesting to check the aggregate stress changes, and how each subjects stress reactivity looks like. For this, we make an average dataframe across the weeks, so we can look at the average change in mood.
#Average Stress per Week:
EMA_avg <- EMA_Data
# Now I can do some merging
EMA_avg <- aggregate(EMA_avg, by=list(EMA_avg$castor_record_id, EMA_avg$Week_Type), FUN=mean, na.rm=T, na.warn=F)
EMA_avg$castor_record_id <- EMA_avg$Group.1
EMA_avg$Week_Type <- EMA_avg$Group.2
EMA_avg$Week_Type <- factor(EMA_avg$Week_Type, levels=c('Control','Stress'),
labels=c('Control','Exam'))
# Split the averaged DF to control and stress week
EMA_avg_exam <- subset(EMA_avg, Week_Type=="Exam")
EMA_avg_control <- subset(EMA_avg, Week_Type!="Exam")
# Merge by week to calculate change
EMA_avg_wide <- merge(EMA_avg_exam, EMA_avg_control, by='castor_record_id', suffixes=c('_exam', '_control'))
# Make Aggregated variables
EMA_avg_wide$mood_positive <-((EMA_avg_wide$mood_positive_control + EMA_avg_wide$mood_positive_exam)/2)
EMA_avg_wide$mood_negative <-((EMA_avg_wide$mood_negative_control + EMA_avg_wide$mood_negative_exam)/2)
EMA_avg_wide$activity_tot <-((EMA_avg_wide$activity_tot_control + EMA_avg_wide$activity_tot_exam)/2)
EMA_avg_wide$social_tot <-((EMA_avg_wide$social_tot_control + EMA_avg_wide$social_tot_exam)/2)
EMA_avg_wide$event_tot <-((EMA_avg_wide$event_tot_control + EMA_avg_wide$event_tot_exam)/2)
EMA_avg_wide$physical_tot <-((EMA_avg_wide$physical_tot_control + EMA_avg_wide$physical_tot_exam)/2)
# Make stress change, and mean centrer it
EMA_avg_wide$stress_reactivity <- EMA_avg_wide$mean_stress_exam - EMA_avg_wide$mean_stress_control
EMA_avg_wide$stress_reactivity_c <- EMA_avg_wide$mean_stress_c_exam - EMA_avg_wide$mean_stress_c_control
EMA_avg_wide$stress_reactivity_z <- scale(EMA_avg_wide$stress_reactivity, center = TRUE, scale = TRUE)
#Box Plot
mean_box <- ggplot(EMA_avg, aes(y=mean_stress, x=Week_Type, colour=Week_Type,fill=Week_Type, na.rm = TRUE)) +
geom_boxplot(alpha=1/2) + geom_jitter(width=0.1, alpha=1/2)+
scale_y_continuous()+ scale_x_discrete()+
ggtitle("Stress Levels per Week")+
xlab("") + ylab("Aggregated stress measure\n")+ ggtheme
mean_box #+ coord_fixed(ratio=1)
We can now run the actual statistics. We use mixed effects models to first validate our paradigm by looking at overall changes in our subjective stress measures between the two weeks. We then go one to investigate the effects of the stress exposure (i.e., stress week) on our outcomes of affect and arousal. We finally investigate and explore the moment-to-moment associations between subjective stress and the outcomes, and then further explore these effects with follow-up analyses and a mediation analysis.
All our models are fit in a maximal approach. We fit fixed effects for all the main effect of interest, with a random slope. We also fit fixed effects with fixed slopes for all the covariates. If a covariate is categorical with fewer than 6 levels, then we only fit fixed effects and no random slopes. Finally, we model subject as the main random effect. For each model, we test the multiple families and links, and using the AIC and residual fit, we select the best model. In the results below, we present the final selected models for brevity. Models are commented when producing this notebook to simplify the generation of the html file.
Additionally, we check models for mulitcolinearity using models without interaction terms. Models reported here are the final ones, with the interaction terms which inflates multicolinearity estimates.
The first set of analysis look at differences between the stress and control weeks. We model this as a fixed effect, and additionally model covariates as fixed effects for sex, study program, week order, day and time, movement, and exercise. We model fixed slopes for the covariates with multiple levels, and a random slope for the fixed effect of interest.
The first models will look at the subjective stress measures, and how those change between the weeks. Based on these models, we can see that there is a significant overall effect of week type on most of our stress measures, which is great. Our stress induction is working. This means we have higher subjective stress in our stress week, indicating our design is achieving its goal.
Control items assessing physical discomfort or illness.
# Model
# glmer.physical_week<- lmer( physical_tot_s ~ Week_Type + #Model week
# Sex + Program + # Model sex and program
# First_Week + ema_day*ema_beep_f + # Model day related items
# acc_delta + physical_excercise_dur + acc_delta*physical_excercise_dur + # Model movement
# (1+Week_Type|castor_record_id) + (0+acc_delta|castor_record_id) +
# (0+ema_survey|castor_record_id) + (0+physical_excercise_dur|castor_record_id) +
# (1|castor_record_id:ema_day:ema_beep_f),
# data=EMA_Data,
# #family=Gamma,
# control=lmerControl(calc.derivs = FALSE))
# Model Summary
model_output(glmer.physical_week)
[4mMODEL INFO:[24m
[3mObservations:[23m 4794
[3mDependent Variable:[23m physical_tot_s
[3mType:[23m Mixed effects linear regression
[4mMODEL FIT:[24m
[3mAIC[23m = 15852.861, [3mBIC[23m = 16221.943
[4mFIXED EFFECTS:
[24m-----------------------------------------------------------------------------------
Est. S.E. t val. d.f. p
-------------------------------------- -------- ------- -------- ---------- -------
(Intercept) 4.593 0.246 18.701 130.233 0.000
Week_TypeStress -0.061 0.094 -0.645 77.815 0.521
SexMale -0.291 0.212 -1.372 77.028 0.174
ProgramBSc_Med -0.078 0.233 -0.336 77.273 0.738
First_WeekExam-First 0.286 0.218 1.312 77.239 0.193
ema_dayEMA 2. -0.225 0.157 -1.433 2717.622 0.152
ema_dayEMA 3. 0.015 0.166 0.092 2745.046 0.927
ema_dayEMA 4. -0.345 0.172 -1.999 2276.572 0.046
ema_dayEMA 5. -0.446 0.177 -2.526 1369.682 0.012
ema_dayEMA 6. -0.280 0.187 -1.498 852.979 0.134
ema_dayEMA 7. -0.440 0.202 -2.184 599.691 0.029
ema_beep_f2 -0.240 0.155 -1.550 2619.050 0.121
ema_beep_f3 -0.328 0.155 -2.117 2695.627 0.034
ema_beep_f4 -0.005 0.154 -0.034 2658.061 0.973
ema_beep_f5 -0.809 0.165 -4.895 2856.021 0.000
ema_beep_f6 -0.612 0.158 -3.873 2779.392 0.000
acc_delta -0.014 0.003 -4.112 4567.600 0.000
physical_excercise_dur 0.607 0.180 3.382 133.755 0.001
ema_dayEMA 2.:ema_beep_f2 0.042 0.214 0.196 2494.454 0.845
ema_dayEMA 3.:ema_beep_f2 -0.081 0.221 -0.364 2577.349 0.716
ema_dayEMA 4.:ema_beep_f2 -0.018 0.224 -0.078 2698.537 0.938
ema_dayEMA 5.:ema_beep_f2 0.153 0.219 0.699 2614.609 0.484
ema_dayEMA 6.:ema_beep_f2 -0.044 0.221 -0.199 2583.679 0.843
ema_dayEMA 7.:ema_beep_f2 0.152 0.228 0.667 2623.559 0.505
ema_dayEMA 2.:ema_beep_f3 0.179 0.216 0.830 2511.717 0.406
ema_dayEMA 3.:ema_beep_f3 -0.264 0.220 -1.202 2605.486 0.230
ema_dayEMA 4.:ema_beep_f3 0.067 0.223 0.302 2714.789 0.763
ema_dayEMA 5.:ema_beep_f3 0.231 0.220 1.051 2637.320 0.293
ema_dayEMA 6.:ema_beep_f3 -0.117 0.221 -0.530 2594.353 0.596
ema_dayEMA 7.:ema_beep_f3 0.067 0.226 0.299 2682.454 0.765
ema_dayEMA 2.:ema_beep_f4 0.055 0.215 0.256 2523.859 0.798
ema_dayEMA 3.:ema_beep_f4 -0.298 0.218 -1.368 2563.772 0.171
ema_dayEMA 4.:ema_beep_f4 -0.043 0.221 -0.194 2654.929 0.846
ema_dayEMA 5.:ema_beep_f4 0.275 0.220 1.250 2643.635 0.211
ema_dayEMA 6.:ema_beep_f4 -0.199 0.223 -0.893 2635.061 0.372
ema_dayEMA 7.:ema_beep_f4 0.269 0.230 1.167 2722.476 0.243
ema_dayEMA 2.:ema_beep_f5 0.172 0.222 0.775 2627.104 0.438
ema_dayEMA 3.:ema_beep_f5 -0.037 0.230 -0.160 2742.037 0.873
ema_dayEMA 4.:ema_beep_f5 0.153 0.229 0.666 2773.748 0.506
ema_dayEMA 5.:ema_beep_f5 0.498 0.236 2.116 2851.437 0.034
ema_dayEMA 6.:ema_beep_f5 0.254 0.236 1.077 2769.162 0.282
ema_dayEMA 7.:ema_beep_f5 0.515 0.237 2.173 2864.959 0.030
ema_dayEMA 2.:ema_beep_f6 0.203 0.219 0.928 2604.480 0.353
ema_dayEMA 3.:ema_beep_f6 -0.200 0.223 -0.894 2673.127 0.372
ema_dayEMA 4.:ema_beep_f6 0.058 0.223 0.260 2741.291 0.795
ema_dayEMA 5.:ema_beep_f6 0.284 0.223 1.269 2717.177 0.205
ema_dayEMA 6.:ema_beep_f6 0.342 0.228 1.503 2686.565 0.133
ema_dayEMA 7.:ema_beep_f6 0.291 0.260 1.123 3030.450 0.262
acc_delta:physical_excercise_dur -0.013 0.024 -0.567 982.327 0.571
-----------------------------------------------------------------------------------
[3m[36mp values calculated using Satterthwaite d.f.[39m[23m
[4m
RANDOM EFFECTS:
[24m--------------------------------------------------------------------------
Group Parameter Std. Dev.
------------------------------------- ------------------------ -----------
castor_record_id.ema_day.ema_beep_f (Intercept) 0.222
castor_record_id physical_excercise_dur 0.246
castor_record_id.1 ema_survey 0.027
castor_record_id.2 acc_delta 0.000
castor_record_id.3 (Intercept) 0.936
castor_record_id.3 Week_TypeStress 0.788
Residual 1.144
--------------------------------------------------------------------------
[4m
Grouping variables:
[24m--------------------------------------------------------
Group # groups ICC
------------------------------------- ---------- -------
castor_record_id:ema_day:ema_beep_f 2962 0.022
castor_record_id 82 0.026
--------------------------------------------------------
[34m# Check for Multicollinearity
[39m
[32mLow Correlation
[39m Parameter VIF Increased SE
Week_Type 1.00 1.00
Sex 1.02 1.01
Program 1.02 1.01
First_Week 1.01 1.00
acc_delta 1.15 1.07
physical_excercise_dur 2.00 1.42
acc_delta:physical_excercise_dur 2.13 1.46
[31mHigh Correlation
[39m Parameter VIF Increased SE
ema_day 21246.19 145.76
ema_beep_f 12290.46 110.86
ema_day:ema_beep_f 81868656.53 9048.13
We tabulate the results, and produce a regression effect plot to visualize the results.
tab_model(glmer.event_week, glmer.activity_week, glmer.social_week, glmer.physical_week,
terms=c("(Intercept)","Week_Type [Stress]", "Sex [Male]", "Program [BSc_Med]", "First_Week [Exam-First]", "acc_delta", "physical_excercise_dur"),
dv.labels=c("Event", "Activitiy", "Social", "Physical"),
title="Table 1. Subjective Stress vs Week Type", show.df=T, show.fstat = T,
transform=NULL,
show.stat=TRUE,
show.se=TRUE) %>%
return() %$%
knitr %>%
asis_output()
| Event | Activitiy | Social | Physical | |||||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | Statistic | p | df | Estimates | std. Error | CI | Statistic | p | df | Estimates | std. Error | CI | Statistic | p | df | Estimates | std. Error | CI | Statistic | p | df |
| (Intercept) | 5.57 | 0.15 | 5.28 – 5.85 | 37.78 | <0.001 | 4737.00 | 4.00 | 0.24 | 3.52 – 4.47 | 16.46 | <0.001 | 4737.00 | 1.46 | 0.05 | 1.37 – 1.55 | 32.05 | <0.001 | 4737.00 | 4.59 | 0.25 | 4.11 – 5.07 | 18.70 | <0.001 | 4737.00 |
| Week_Type [Stress] | 0.31 | 0.06 | 0.18 – 0.43 | 4.90 | <0.001 | 4737.00 | 0.49 | 0.10 | 0.29 – 0.69 | 4.81 | <0.001 | 4737.00 | 0.02 | 0.02 | -0.01 – 0.05 | 1.19 | 0.236 | 4737.00 | -0.06 | 0.09 | -0.25 – 0.12 | -0.65 | 0.519 | 4737.00 |
| Sex [Male] | -0.01 | 0.09 | -0.19 – 0.17 | -0.11 | 0.910 | 4737.00 | -0.05 | 0.16 | -0.36 – 0.26 | -0.32 | 0.752 | 4737.00 | 0.08 | 0.02 | 0.03 – 0.12 | 3.09 | 0.002 | 4737.00 | -0.29 | 0.21 | -0.71 – 0.12 | -1.37 | 0.170 | 4737.00 |
| Program [BSc_Med] | -0.01 | 0.10 | -0.21 – 0.18 | -0.15 | 0.881 | 4737.00 | 0.00 | 0.17 | -0.33 – 0.34 | 0.02 | 0.980 | 4737.00 | 0.01 | 0.03 | -0.04 – 0.06 | 0.29 | 0.770 | 4737.00 | -0.08 | 0.23 | -0.54 – 0.38 | -0.34 | 0.737 | 4737.00 |
| First_Week [Exam-First] | -0.07 | 0.09 | -0.25 – 0.11 | -0.79 | 0.431 | 4737.00 | 0.10 | 0.16 | -0.22 – 0.42 | 0.62 | 0.538 | 4737.00 | 0.03 | 0.03 | -0.02 – 0.08 | 1.08 | 0.282 | 4737.00 | 0.29 | 0.22 | -0.14 – 0.71 | 1.31 | 0.190 | 4737.00 |
| acc_delta | -0.00 | 0.00 | -0.01 – 0.01 | -0.26 | 0.793 | 4737.00 | 0.00 | 0.01 | -0.01 – 0.01 | 0.15 | 0.881 | 4737.00 | -0.00 | 0.00 | -0.00 – 0.00 | -1.34 | 0.179 | 4737.00 | -0.01 | 0.00 | -0.02 – -0.01 | -4.11 | <0.001 | 4737.00 |
| physical_excercise_dur | -0.30 | 0.18 | -0.66 – 0.06 | -1.63 | 0.104 | 4737.00 | 0.17 | 0.27 | -0.36 – 0.70 | 0.63 | 0.532 | 4737.00 | -0.08 | 0.06 | -0.20 – 0.04 | -1.36 | 0.174 | 4737.00 | 0.61 | 0.18 | 0.26 – 0.96 | 3.38 | 0.001 | 4737.00 |
| Random Effects | ||||||||||||||||||||||||
| σ2 | 1.31 | |||||||||||||||||||||||
| τ00 | 0.05 castor_record_id.ema_day.ema_beep_f | |||||||||||||||||||||||
| 0.06 castor_record_id | ||||||||||||||||||||||||
| 0.00 castor_record_id.1 | ||||||||||||||||||||||||
| 0.00 castor_record_id.2 | ||||||||||||||||||||||||
| 0.88 castor_record_id.3 | ||||||||||||||||||||||||
| τ11 | 0.62 castor_record_id.3.Week_TypeStress | |||||||||||||||||||||||
| ρ01 | -0.35 castor_record_id.3 | |||||||||||||||||||||||
| N | 82 castor_record_id | 82 castor_record_id | 82 castor_record_id | 82 castor_record_id | ||||||||||||||||||||
| 7 ema_day | 7 ema_day | 7 ema_day | 7 ema_day | |||||||||||||||||||||
| 6 ema_beep_f | 6 ema_beep_f | 6 ema_beep_f | 6 ema_beep_f | |||||||||||||||||||||
| Observations | 4794 | 4794 | 4794 | 4794 | ||||||||||||||||||||
| Marginal R2 / Conditional R2 | NA | NA | NA | 0.083 / NA | ||||||||||||||||||||
# Subset the variables we want to plot
plot.mood_week <- plot_models(glmer.event_week, glmer.activity_week, glmer.social_week, glmer.physical_week)
rm_term <- as.vector(plot.mood_week$data$term[1:length(plot.mood_week$data$term)])
rm_term <- rm_term[rm_term != "Week_TypeStress"]
# Make the plot
plot.mood_week <- plot_models(glmer.event_week, glmer.activity_week, glmer.social_week, glmer.physical_week,
m.labels=c("Event Stress", "Activity Stress","Social Stress", "Physical Stress" ),
axis.labels = c(" "),
# Stastistical Stuff
rm.terms = rm_term,
#show.values = T,
#value.size = 4,
#std.est=T,
show.p=T,
p.shape=T,
legend.pval.title = "Significance",
# Visual Stuff
colors="#666666",
dot.size=2,
line.size = 1,
spacing=1 ,
vline.color = "darkgrey",
legend.title = "") + ylab("\nParameter Estimate (a.u.)") + xlab("Subjective Stress") +
ptheme +
theme(axis.ticks.y=element_blank(),
legend.text=element_text(size=16), legend.title = element_text(size=16),
axis.title = element_text(size=16),axis.text.x=element_text(size=11),
panel.grid.major.x = element_line(colour = "grey95"),
panel.grid.minor.x = element_line(colour = "grey90"))
# Show and save
plot.mood_week+ coord_flip(ylim=c(-0.65,0.65)) + theme(panel.grid.major.y = element_line(colour = "grey95"), panel.grid.minor.y = element_line(colour = "grey90"))
Coordinate system already present. Adding new coordinate system, which will replace the existing one.
ggsave("figures/fig_WeekResid_SubjectiveStress.tiff",units="in", width=4, height=4, dpi=300, compression = 'lzw', bg="transparent")
TIFFOpen: figures/fig_WeekResid_SubjectiveStress.tiff: No such file or directory.
Next, we want to investigate the impact of an examination period on mood. Conceptually, mood is an outcome of being stressed, and not a direct measure of stress. Below are the results of the models for both positive and negative mood.
# Model
# glmer.posmood_week <- glmer( mood_positive_s ~ Week_Type + #Model week
# Sex + Program + # Model sex and program
# First_Week + ema_day*ema_beep_f + # Model day related items
# physical_excercise_dur + acc_delta + acc_delta*physical_excercise_dur + # Model movement
# (1+Week_Type|castor_record_id) + (0+acc_delta|castor_record_id) +
# (0+ema_survey|castor_record_id) + (0+physical_excercise_dur|castor_record_id) +
# (1|castor_record_id:ema_day:ema_beep_f),
# data=EMA_Data,
# family=gaussian(link="log"),
# control=glmerControl(calc.derivs = FALSE))
# Model Summary
model_output(glmer.posmood_week)
[4mMODEL INFO:[24m
[3mObservations:[23m 4794
[3mDependent Variable:[23m mood_positive_s
[3mType:[23m Mixed effects generalized linear regression
[3mError Distribution: [23mgaussian
[3mLink function: [23mlog
[4mMODEL FIT:[24m
[3mAIC[23m = 16233.449, [3mBIC[23m = 16602.531
[4mFIXED EFFECTS:
[24m------------------------------------------------------------------------
Est. S.E. t val. p
-------------------------------------- -------- ------- -------- -------
(Intercept) 1.825 0.045 40.447 0.000
Week_TypeStress -0.078 0.015 -5.232 0.000
SexMale 0.046 0.041 1.137 0.256
ProgramBSc_Med 0.024 0.045 0.530 0.596
First_WeekExam-First -0.008 0.042 -0.180 0.857
ema_dayEMA 2. 0.036 0.025 1.433 0.152
ema_dayEMA 3. -0.016 0.027 -0.596 0.551
ema_dayEMA 4. -0.011 0.027 -0.401 0.689
ema_dayEMA 5. -0.004 0.027 -0.142 0.887
ema_dayEMA 6. -0.018 0.028 -0.644 0.519
ema_dayEMA 7. -0.029 0.030 -0.980 0.327
ema_beep_f2 0.016 0.025 0.622 0.534
ema_beep_f3 0.031 0.025 1.250 0.211
ema_beep_f4 0.020 0.025 0.808 0.419
ema_beep_f5 0.084 0.026 3.258 0.001
ema_beep_f6 0.076 0.025 3.060 0.002
physical_excercise_dur 0.063 0.027 2.325 0.020
acc_delta 0.001 0.001 2.050 0.040
ema_dayEMA 2.:ema_beep_f2 -0.024 0.034 -0.711 0.477
ema_dayEMA 3.:ema_beep_f2 0.014 0.036 0.381 0.703
ema_dayEMA 4.:ema_beep_f2 0.018 0.036 0.491 0.623
ema_dayEMA 5.:ema_beep_f2 -0.002 0.036 -0.066 0.947
ema_dayEMA 6.:ema_beep_f2 0.000 0.036 0.012 0.991
ema_dayEMA 7.:ema_beep_f2 0.005 0.037 0.147 0.884
ema_dayEMA 2.:ema_beep_f3 -0.075 0.034 -2.172 0.030
ema_dayEMA 3.:ema_beep_f3 -0.016 0.035 -0.463 0.643
ema_dayEMA 4.:ema_beep_f3 0.003 0.036 0.095 0.924
ema_dayEMA 5.:ema_beep_f3 -0.020 0.036 -0.565 0.572
ema_dayEMA 6.:ema_beep_f3 -0.012 0.036 -0.326 0.745
ema_dayEMA 7.:ema_beep_f3 -0.048 0.037 -1.289 0.198
ema_dayEMA 2.:ema_beep_f4 -0.031 0.034 -0.904 0.366
ema_dayEMA 3.:ema_beep_f4 0.038 0.035 1.086 0.277
ema_dayEMA 4.:ema_beep_f4 0.021 0.036 0.583 0.560
ema_dayEMA 5.:ema_beep_f4 -0.019 0.036 -0.536 0.592
ema_dayEMA 6.:ema_beep_f4 -0.007 0.036 -0.197 0.844
ema_dayEMA 7.:ema_beep_f4 -0.003 0.038 -0.092 0.927
ema_dayEMA 2.:ema_beep_f5 -0.050 0.035 -1.452 0.147
ema_dayEMA 3.:ema_beep_f5 -0.009 0.036 -0.235 0.814
ema_dayEMA 4.:ema_beep_f5 0.000 0.036 0.011 0.991
ema_dayEMA 5.:ema_beep_f5 -0.072 0.038 -1.927 0.054
ema_dayEMA 6.:ema_beep_f5 -0.055 0.038 -1.471 0.141
ema_dayEMA 7.:ema_beep_f5 -0.033 0.038 -0.863 0.388
ema_dayEMA 2.:ema_beep_f6 -0.025 0.034 -0.734 0.463
ema_dayEMA 3.:ema_beep_f6 -0.029 0.035 -0.825 0.409
ema_dayEMA 4.:ema_beep_f6 -0.008 0.035 -0.229 0.819
ema_dayEMA 5.:ema_beep_f6 -0.022 0.035 -0.629 0.529
ema_dayEMA 6.:ema_beep_f6 -0.046 0.036 -1.265 0.206
ema_dayEMA 7.:ema_beep_f6 -0.010 0.041 -0.231 0.817
physical_excercise_dur:acc_delta -0.002 0.004 -0.435 0.664
------------------------------------------------------------------------
[4m
RANDOM EFFECTS:
[24m--------------------------------------------------------------------------
Group Parameter Std. Dev.
------------------------------------- ------------------------ -----------
castor_record_id.ema_day.ema_beep_f (Intercept) 0.000
castor_record_id physical_excercise_dur 0.000
castor_record_id.1 ema_survey 0.002
castor_record_id.2 acc_delta 0.000
castor_record_id.3 (Intercept) 0.173
castor_record_id.3 Week_TypeStress 0.122
Residual 1.220
--------------------------------------------------------------------------
[4m
Grouping variables:
[24m--------------------------------------------------------
Group # groups ICC
------------------------------------- ---------- -------
castor_record_id:ema_day:ema_beep_f 2962 0.000
castor_record_id 82 0.000
--------------------------------------------------------
[34m# Check for Multicollinearity
[39m
[32mLow Correlation
[39m Parameter VIF Increased SE
Week_Type 1.00 1.00
Sex 1.02 1.01
Program 1.02 1.01
First_Week 1.01 1.00
physical_excercise_dur 2.09 1.45
acc_delta 1.15 1.07
physical_excercise_dur:acc_delta 2.22 1.49
[31mHigh Correlation
[39m Parameter VIF Increased SE
ema_day 44880.88 211.85
ema_beep_f 10326.16 101.62
ema_day:ema_beep_f 130049717.05 11403.93
# Model
# glmer.negmood_week <- glmer( mood_negative_s ~ Week_Type + #Model week
# Sex + Program + # Model sex and program
# First_Week + ema_day*ema_beep_f + # Model day related items
# physical_excercise_dur + acc_delta + acc_delta*physical_excercise_dur + # Model movement
# (1+Week_Type|castor_record_id) + (0+acc_delta|castor_record_id) +
# (0+ema_survey|castor_record_id) + (0+physical_excercise_dur|castor_record_id) +
# (1|castor_record_id:ema_day:ema_beep_f),
# data=EMA_Data,
# family=Gamma(link="log"),
# control=glmerControl(calc.derivs = FALSE))
# Model Summary
model_output(glmer.negmood_week)
[4mMODEL INFO:[24m
[3mObservations:[23m 4794
[3mDependent Variable:[23m mood_negative_s
[3mType:[23m Mixed effects generalized linear regression
[3mError Distribution: [23mGamma
[3mLink function: [23mlog
[4mMODEL FIT:[24m
[3mAIC[23m = 11278.717, [3mBIC[23m = 11647.799
[4mFIXED EFFECTS:
[24m------------------------------------------------------------------------
Est. S.E. t val. p
-------------------------------------- -------- ------- -------- -------
(Intercept) 1.009 0.070 14.325 0.000
Week_TypeStress 0.113 0.023 4.799 0.000
SexMale -0.188 0.054 -3.500 0.000
ProgramBSc_Med -0.135 0.059 -2.288 0.022
First_WeekExam-First 0.057 0.055 1.027 0.305
ema_dayEMA 2. -0.184 0.060 -3.080 0.002
ema_dayEMA 3. -0.152 0.062 -2.452 0.014
ema_dayEMA 4. -0.183 0.062 -2.936 0.003
ema_dayEMA 5. -0.165 0.062 -2.663 0.008
ema_dayEMA 6. -0.148 0.064 -2.326 0.020
ema_dayEMA 7. -0.107 0.066 -1.629 0.103
ema_beep_f2 -0.018 0.059 -0.312 0.755
ema_beep_f3 -0.132 0.059 -2.243 0.025
ema_beep_f4 -0.136 0.059 -2.306 0.021
ema_beep_f5 -0.221 0.062 -3.543 0.000
ema_beep_f6 -0.287 0.060 -4.790 0.000
physical_excercise_dur -0.055 0.074 -0.746 0.456
acc_delta -0.003 0.001 -2.319 0.020
ema_dayEMA 2.:ema_beep_f2 0.059 0.083 0.717 0.473
ema_dayEMA 3.:ema_beep_f2 0.023 0.085 0.274 0.784
ema_dayEMA 4.:ema_beep_f2 0.006 0.086 0.069 0.945
ema_dayEMA 5.:ema_beep_f2 0.010 0.084 0.115 0.909
ema_dayEMA 6.:ema_beep_f2 -0.007 0.085 -0.086 0.931
ema_dayEMA 7.:ema_beep_f2 -0.060 0.087 -0.682 0.495
ema_dayEMA 2.:ema_beep_f3 0.185 0.083 2.226 0.026
ema_dayEMA 3.:ema_beep_f3 0.110 0.084 1.299 0.194
ema_dayEMA 4.:ema_beep_f3 0.119 0.085 1.398 0.162
ema_dayEMA 5.:ema_beep_f3 0.118 0.084 1.400 0.161
ema_dayEMA 6.:ema_beep_f3 0.106 0.085 1.247 0.212
ema_dayEMA 7.:ema_beep_f3 0.171 0.086 1.989 0.047
ema_dayEMA 2.:ema_beep_f4 0.196 0.083 2.355 0.019
ema_dayEMA 3.:ema_beep_f4 0.091 0.084 1.084 0.279
ema_dayEMA 4.:ema_beep_f4 0.166 0.085 1.963 0.050
ema_dayEMA 5.:ema_beep_f4 0.155 0.084 1.839 0.066
ema_dayEMA 6.:ema_beep_f4 0.142 0.085 1.667 0.095
ema_dayEMA 7.:ema_beep_f4 0.137 0.088 1.566 0.117
ema_dayEMA 2.:ema_beep_f5 0.200 0.085 2.353 0.019
ema_dayEMA 3.:ema_beep_f5 0.092 0.088 1.052 0.293
ema_dayEMA 4.:ema_beep_f5 0.229 0.087 2.626 0.009
ema_dayEMA 5.:ema_beep_f5 0.214 0.089 2.407 0.016
ema_dayEMA 6.:ema_beep_f5 0.185 0.090 2.066 0.039
ema_dayEMA 7.:ema_beep_f5 0.161 0.090 1.799 0.072
ema_dayEMA 2.:ema_beep_f6 0.258 0.084 3.083 0.002
ema_dayEMA 3.:ema_beep_f6 0.245 0.085 2.866 0.004
ema_dayEMA 4.:ema_beep_f6 0.232 0.085 2.732 0.006
ema_dayEMA 5.:ema_beep_f6 0.282 0.085 3.314 0.001
ema_dayEMA 6.:ema_beep_f6 0.303 0.087 3.489 0.000
ema_dayEMA 7.:ema_beep_f6 0.314 0.097 3.222 0.001
physical_excercise_dur:acc_delta -0.004 0.008 -0.529 0.597
------------------------------------------------------------------------
[4m
RANDOM EFFECTS:
[24m--------------------------------------------------------------------------
Group Parameter Std. Dev.
------------------------------------- ------------------------ -----------
castor_record_id.ema_day.ema_beep_f (Intercept) 0.250
castor_record_id physical_excercise_dur 0.338
castor_record_id.1 ema_survey 0.004
castor_record_id.2 acc_delta 0.006
castor_record_id.3 (Intercept) 0.241
castor_record_id.3 Week_TypeStress 0.190
Residual 0.315
--------------------------------------------------------------------------
[4m
Grouping variables:
[24m--------------------------------------------------------
Group # groups ICC
------------------------------------- ---------- -------
castor_record_id:ema_day:ema_beep_f 2962 0.188
castor_record_id 82 0.342
--------------------------------------------------------
[34m# Check for Multicollinearity
[39m
[32mLow Correlation
[39m Parameter VIF Increased SE
Week_Type 1.00 1.00
Sex 1.02 1.01
Program 1.02 1.01
First_Week 1.01 1.00
physical_excercise_dur 1.52 1.23
acc_delta 1.12 1.06
physical_excercise_dur:acc_delta 1.63 1.28
[31mHigh Correlation
[39m Parameter VIF Increased SE
ema_day 39883.11 199.71
ema_beep_f 12266.18 110.75
ema_day:ema_beep_f 140648388.77 11859.53
tab_model(glmer.posmood_week, glmer.negmood_week,
terms=c("(Intercept)","Week_Type [Stress]", "Sex [Male]", "Program [BSc_Med]", "First_Week [Exam-First]", "acc_delta", "physical_excercise_dur"),
dv.labels=c("Positive", "Negative"),
title="Table 2. Mood vs Week Type",
transform=NULL,
show.stat=TRUE,
show.se=TRUE) %>%
return() %$%
knitr %>%
asis_output()
| Positive | Negative | |||||||||
|---|---|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | Statistic | p | Estimates | std. Error | CI | Statistic | p |
| (Intercept) | 1.83 | 0.05 | 1.74 – 1.91 | 40.45 | <0.001 | 1.01 | 0.07 | 0.87 – 1.15 | 14.32 | <0.001 |
| Week_Type [Stress] | -0.08 | 0.01 | -0.11 – -0.05 | -5.23 | <0.001 | 0.11 | 0.02 | 0.07 – 0.16 | 4.80 | <0.001 |
| Sex [Male] | 0.05 | 0.04 | -0.03 – 0.13 | 1.14 | 0.256 | -0.19 | 0.05 | -0.29 – -0.08 | -3.50 | <0.001 |
| Program [BSc_Med] | 0.02 | 0.04 | -0.06 – 0.11 | 0.53 | 0.596 | -0.14 | 0.06 | -0.25 – -0.02 | -2.29 | 0.022 |
| First_Week [Exam-First] | -0.01 | 0.04 | -0.09 – 0.07 | -0.18 | 0.857 | 0.06 | 0.06 | -0.05 – 0.16 | 1.03 | 0.305 |
| physical_excercise_dur | 0.06 | 0.03 | 0.01 – 0.12 | 2.32 | 0.020 | -0.06 | 0.07 | -0.20 – 0.09 | -0.75 | 0.456 |
| acc_delta | 0.00 | 0.00 | 0.00 – 0.00 | 2.05 | 0.040 | -0.00 | 0.00 | -0.01 – -0.00 | -2.32 | 0.020 |
| Random Effects | ||||||||||
| σ2 | 1.49 | |||||||||
| τ00 | 0.00 castor_record_id.ema_day.ema_beep_f | |||||||||
| 0.00 castor_record_id | ||||||||||
| 0.00 castor_record_id.1 | ||||||||||
| 0.00 castor_record_id.2 | ||||||||||
| 0.03 castor_record_id.3 | ||||||||||
| τ11 | 0.01 castor_record_id.3.Week_TypeStress | |||||||||
| ρ01 | -0.15 castor_record_id.3 | |||||||||
| N | 82 castor_record_id | 82 castor_record_id | ||||||||
| 7 ema_day | 7 ema_day | |||||||||
| 6 ema_beep_f | 6 ema_beep_f | |||||||||
| Observations | 4794 | 4794 | ||||||||
| Marginal R2 / Conditional R2 | 0.002 / NA | NA | ||||||||
We established changes in subjective stress between the weeks, which has an impact on our mood outcomes. We also want to see if there’s a general difference between the physiology variables and the weeks.To do this we apply the same approach for the physiology measures as the other models. We also add measures of skin temperature here however, to make sure that we factor this in. Skin temperature is important as it can be related to sweat or heat, which can also lead to increased arousal measures.
First we model the skin conductance, controlling for skin temperature and physical exercise.
# Model Summary
model_output(glmer.sc_ton_mean_week)
Could not calculate r-squared. Try removing missing data
before fitting the model.
[4mMODEL INFO:[24m
[3mObservations:[23m 4565
[3mDependent Variable:[23m sc_tonic_mean_s
[3mType:[23m Mixed effects generalized linear regression
[3mError Distribution: [23mGamma
[3mLink function: [23mlog
[4mMODEL FIT:[24m
[3mAIC[23m = -1329.642, [3mBIC[23m = -950.497
[4mFIXED EFFECTS:
[24m-----------------------------------------------------------------
Est. S.E. t val. p
------------------------------- -------- ------- -------- -------
(Intercept) 0.086 0.031 2.786 0.005
Week_TypeStress -0.004 0.011 -0.335 0.738
SexMale -0.007 0.013 -0.552 0.581
ProgramBSc_Med 0.033 0.015 2.256 0.024
First_WeekExam-First -0.038 0.014 -2.762 0.006
ema_beep_f2 -0.027 0.037 -0.740 0.460
ema_beep_f3 -0.021 0.036 -0.580 0.562
ema_beep_f4 -0.045 0.036 -1.242 0.214
ema_beep_f5 -0.027 0.038 -0.716 0.474
ema_beep_f6 0.011 0.037 0.292 0.771
ema_dayEMA 2. -0.026 0.037 -0.701 0.483
ema_dayEMA 3. -0.077 0.038 -2.025 0.043
ema_dayEMA 4. -0.035 0.038 -0.912 0.362
ema_dayEMA 5. -0.097 0.038 -2.573 0.010
ema_dayEMA 6. -0.068 0.039 -1.763 0.078
ema_dayEMA 7. -0.063 0.039 -1.621 0.105
physical_excercise_dur 0.076 0.037 2.071 0.038
acc_delta 0.014 0.001 10.513 0.000
temp_mean_z 0.017 0.007 2.379 0.017
temp_slope_z 0.019 0.008 2.342 0.019
ema_beep_f2:ema_dayEMA 2. -0.002 0.051 -0.043 0.966
ema_beep_f3:ema_dayEMA 2. 0.039 0.051 0.757 0.449
ema_beep_f4:ema_dayEMA 2. 0.027 0.051 0.527 0.598
ema_beep_f5:ema_dayEMA 2. -0.030 0.052 -0.570 0.569
ema_beep_f6:ema_dayEMA 2. -0.031 0.052 -0.591 0.555
ema_beep_f2:ema_dayEMA 3. 0.026 0.052 0.497 0.619
ema_beep_f3:ema_dayEMA 3. 0.022 0.052 0.431 0.666
ema_beep_f4:ema_dayEMA 3. 0.052 0.052 1.006 0.315
ema_beep_f5:ema_dayEMA 3. 0.063 0.054 1.170 0.242
ema_beep_f6:ema_dayEMA 3. 0.058 0.052 1.106 0.269
ema_beep_f2:ema_dayEMA 4. -0.014 0.053 -0.260 0.795
ema_beep_f3:ema_dayEMA 4. -0.039 0.052 -0.743 0.457
ema_beep_f4:ema_dayEMA 4. -0.011 0.052 -0.205 0.838
ema_beep_f5:ema_dayEMA 4. -0.008 0.053 -0.158 0.875
ema_beep_f6:ema_dayEMA 4. -0.044 0.052 -0.848 0.397
ema_beep_f2:ema_dayEMA 5. 0.063 0.052 1.204 0.229
ema_beep_f3:ema_dayEMA 5. 0.046 0.052 0.878 0.380
ema_beep_f4:ema_dayEMA 5. 0.060 0.052 1.151 0.250
ema_beep_f5:ema_dayEMA 5. 0.035 0.055 0.633 0.527
ema_beep_f6:ema_dayEMA 5. 0.042 0.052 0.803 0.422
ema_beep_f2:ema_dayEMA 6. 0.047 0.053 0.889 0.374
ema_beep_f3:ema_dayEMA 6. 0.016 0.052 0.306 0.760
ema_beep_f4:ema_dayEMA 6. 0.026 0.053 0.501 0.617
ema_beep_f5:ema_dayEMA 6. 0.019 0.055 0.341 0.733
ema_beep_f6:ema_dayEMA 6. -0.008 0.054 -0.151 0.880
ema_beep_f2:ema_dayEMA 7. 0.026 0.054 0.480 0.631
ema_beep_f3:ema_dayEMA 7. 0.038 0.053 0.719 0.472
ema_beep_f4:ema_dayEMA 7. 0.071 0.054 1.311 0.190
ema_beep_f5:ema_dayEMA 7. 0.040 0.055 0.725 0.469
ema_beep_f6:ema_dayEMA 7. -0.024 0.059 -0.406 0.685
-----------------------------------------------------------------
[4m
RANDOM EFFECTS:
[24m--------------------------------------------------------------------------
Group Parameter Std. Dev.
------------------------------------- ------------------------ -----------
castor_record_id.ema_day.ema_beep_f (Intercept) 0.160
castor_record_id acc_delta 0.010
castor_record_id.1 physical_excercise_dur 0.208
castor_record_id.2 temp_mean_z 0.046
castor_record_id.3 temp_slope_z 0.056
castor_record_id.4 (Intercept) 0.058
castor_record_id.4 Week_TypeStress 0.077
Residual 0.175
--------------------------------------------------------------------------
[4m
Grouping variables:
[24m--------------------------------------------------------
Group # groups ICC
------------------------------------- ---------- -------
castor_record_id:ema_day:ema_beep_f 2879 0.238
castor_record_id 82 0.001
--------------------------------------------------------
Model has interaction terms. VIFs might be inflated. You may check multicollinearity among predictors of a model without interaction terms.
[34m# Check for Multicollinearity
[39m
[32mLow Correlation
[39m Parameter VIF Increased SE
Week_Type 1.00 1.00
Sex 1.02 1.01
Program 1.02 1.01
First_Week 1.00 1.00
physical_excercise_dur 1.01 1.00
acc_delta 1.01 1.01
temp_mean_z 1.02 1.01
temp_slope_z 1.01 1.00
[31mHigh Correlation
[39m Parameter VIF Increased SE
ema_beep_f 11831.31 108.77
ema_day 55890.06 236.41
ema_beep_f:ema_day 181695716.16 13479.46
# Model
# glmer.sc_ph_num <- glmer( sc_phasic_num ~ Week_Type +
# Sex + Program + # Model pop differences
# First_Week + ema_beep_f*ema_day + # Model day related differencxes
# physical_excercise_dur + acc_delta + # Modle movement stuff
# temp_mean_z + temp_slope_z + # Model temp effects
# (1+Week_Type |castor_record_id) +
# (0+temp_slope_z|castor_record_id) + (0+temp_mean_z|castor_record_id) +
# (0+physical_excercise_dur|castor_record_id) + (0+acc_delta|castor_record_id) +
# (1 | castor_record_id:ema_day:ema_beep_f),
# EMA_Data,
# family=poisson(link="log"),
# control=lmerControl(calc.derivs = FALSE) )
# Model Summary
model_output(glmer.sc_ph_num)
[4mMODEL INFO:[24m
[3mObservations:[23m 4565
[3mDependent Variable:[23m sc_phasic_num
[3mType:[23m Mixed effects generalized linear regression
[3mError Distribution: [23mpoisson
[3mLink function: [23mlog
[4mMODEL FIT:[24m
[3mAIC[23m = 44355.736, [3mBIC[23m = 44728.454
[4mFIXED EFFECTS:
[24m-----------------------------------------------------------------
Est. S.E. z val. p
------------------------------- -------- ------- -------- -------
(Intercept) 1.368 0.309 4.424 0.000
Week_TypeStress -0.233 0.078 -3.000 0.003
SexMale 0.466 0.257 1.817 0.069
ProgramBSc_Med 0.199 0.283 0.701 0.483
First_WeekExam-First -0.001 0.264 -0.006 0.996
ema_beep_f2 -0.049 0.224 -0.221 0.825
ema_beep_f3 -0.227 0.221 -1.029 0.303
ema_beep_f4 -0.510 0.223 -2.286 0.022
ema_beep_f5 -0.393 0.232 -1.690 0.091
ema_beep_f6 -0.236 0.227 -1.040 0.298
ema_dayEMA 2. -0.028 0.225 -0.127 0.899
ema_dayEMA 3. -0.427 0.231 -1.849 0.065
ema_dayEMA 4. -0.510 0.229 -2.230 0.026
ema_dayEMA 5. -0.564 0.230 -2.451 0.014
ema_dayEMA 6. -0.434 0.236 -1.834 0.067
ema_dayEMA 7. -0.565 0.237 -2.387 0.017
physical_excercise_dur 0.108 0.289 0.374 0.708
acc_delta 0.150 0.014 10.612 0.000
temp_mean_z 0.323 0.067 4.843 0.000
temp_slope_z 0.153 0.058 2.661 0.008
ema_beep_f2:ema_dayEMA 2. -0.381 0.314 -1.211 0.226
ema_beep_f3:ema_dayEMA 2. 0.049 0.314 0.157 0.875
ema_beep_f4:ema_dayEMA 2. 0.273 0.316 0.864 0.387
ema_beep_f5:ema_dayEMA 2. -0.042 0.322 -0.129 0.897
ema_beep_f6:ema_dayEMA 2. -0.640 0.322 -1.987 0.047
ema_beep_f2:ema_dayEMA 3. 0.180 0.321 0.560 0.576
ema_beep_f3:ema_dayEMA 3. 0.282 0.318 0.887 0.375
ema_beep_f4:ema_dayEMA 3. 0.480 0.317 1.514 0.130
ema_beep_f5:ema_dayEMA 3. 0.326 0.329 0.990 0.322
ema_beep_f6:ema_dayEMA 3. 0.235 0.323 0.726 0.468
ema_beep_f2:ema_dayEMA 4. 0.126 0.322 0.390 0.697
ema_beep_f3:ema_dayEMA 4. 0.172 0.320 0.536 0.592
ema_beep_f4:ema_dayEMA 4. 0.573 0.320 1.791 0.073
ema_beep_f5:ema_dayEMA 4. 0.233 0.327 0.714 0.475
ema_beep_f6:ema_dayEMA 4. 0.109 0.322 0.339 0.735
ema_beep_f2:ema_dayEMA 5. -0.139 0.321 -0.433 0.665
ema_beep_f3:ema_dayEMA 5. 0.185 0.320 0.580 0.562
ema_beep_f4:ema_dayEMA 5. 0.673 0.319 2.107 0.035
ema_beep_f5:ema_dayEMA 5. -0.011 0.339 -0.033 0.973
ema_beep_f6:ema_dayEMA 5. -0.270 0.328 -0.824 0.410
ema_beep_f2:ema_dayEMA 6. 0.017 0.325 0.053 0.958
ema_beep_f3:ema_dayEMA 6. 0.256 0.324 0.790 0.430
ema_beep_f4:ema_dayEMA 6. 0.436 0.327 1.334 0.182
ema_beep_f5:ema_dayEMA 6. 0.126 0.340 0.371 0.711
ema_beep_f6:ema_dayEMA 6. -0.209 0.334 -0.625 0.532
ema_beep_f2:ema_dayEMA 7. 0.366 0.327 1.118 0.264
ema_beep_f3:ema_dayEMA 7. 0.351 0.324 1.084 0.278
ema_beep_f4:ema_dayEMA 7. 0.348 0.332 1.050 0.294
ema_beep_f5:ema_dayEMA 7. 0.073 0.336 0.217 0.828
ema_beep_f6:ema_dayEMA 7. 0.099 0.364 0.273 0.785
-----------------------------------------------------------------
[4m
RANDOM EFFECTS:
[24m--------------------------------------------------------------------------
Group Parameter Std. Dev.
------------------------------------- ------------------------ -----------
castor_record_id.ema_day.ema_beep_f (Intercept) 1.262
castor_record_id acc_delta 0.124
castor_record_id.1 physical_excercise_dur 2.320
castor_record_id.2 temp_mean_z 0.575
castor_record_id.3 temp_slope_z 0.489
castor_record_id.4 (Intercept) 1.103
castor_record_id.4 Week_TypeStress 0.678
--------------------------------------------------------------------------
[4m
Grouping variables:
[24m--------------------------------------------------------
Group # groups ICC
------------------------------------- ---------- -------
castor_record_id:ema_day:ema_beep_f 2879 0.163
castor_record_id 82 0.002
--------------------------------------------------------
[34m# Check for Multicollinearity
[39m
[32mLow Correlation
[39m Parameter VIF Increased SE
Week_Type 1.00 1.00
Sex 1.02 1.01
Program 1.02 1.01
First_Week 1.00 1.00
physical_excercise_dur 1.00 1.00
acc_delta 1.00 1.00
temp_mean_z 1.00 1.00
temp_slope_z 1.00 1.00
[31mHigh Correlation
[39m Parameter VIF Increased SE
ema_beep_f 10999.77 104.88
ema_day 46175.74 214.89
ema_beep_f:ema_day 144610265.17 12025.40
# Model
# glmer.sc_ph_mag <- glmer( sc_phasic_mag_s ~ Week_Type +
# Sex + Program + # Model pop differences
# First_Week + ema_beep_f*ema_day + # Model day related differencxes
# physical_excercise_dur +acc_delta + # Modle movement stuff
# temp_mean_z + temp_slope_z + # Model temp effects
# (1+Week_Type |castor_record_id) +
# (0+temp_slope_z|castor_record_id) + (0+temp_mean_z|castor_record_id) +
# (0+physical_excercise_dur|castor_record_id) + (0+acc_delta|castor_record_id) +
# (1 | castor_record_id:ema_day:ema_beep_f),
# EMA_Data,
# family=Gamma(link="identity"),
# control=glmerControl(calc.derivs = FALSE) )
# Model Summary
model_output(glmer.sc_ph_mag)
[4mMODEL INFO:[24m
[3mObservations:[23m 4565
[3mDependent Variable:[23m sc_phasic_mag_s
[3mType:[23m Mixed effects generalized linear regression
[3mError Distribution: [23mGamma
[3mLink function: [23midentity
[4mMODEL FIT:[24m
[3mAIC[23m = 10389.643, [3mBIC[23m = 10768.787
[4mFIXED EFFECTS:
[24m-----------------------------------------------------------------
Est. S.E. t val. p
------------------------------- -------- ------- -------- -------
(Intercept) 2.136 0.122 17.572 0.000
Week_TypeStress -0.103 0.041 -2.493 0.013
SexMale 0.055 0.061 0.892 0.372
ProgramBSc_Med 0.046 0.067 0.685 0.493
First_WeekExam-First -0.013 0.063 -0.203 0.839
ema_beep_f2 -0.103 0.136 -0.759 0.448
ema_beep_f3 -0.130 0.136 -0.956 0.339
ema_beep_f4 -0.168 0.135 -1.246 0.213
ema_beep_f5 -0.198 0.140 -1.411 0.158
ema_beep_f6 -0.027 0.137 -0.195 0.846
ema_dayEMA 2. -0.107 0.137 -0.781 0.435
ema_dayEMA 3. -0.233 0.140 -1.666 0.096
ema_dayEMA 4. -0.126 0.141 -0.894 0.372
ema_dayEMA 5. -0.244 0.137 -1.780 0.075
ema_dayEMA 6. -0.223 0.141 -1.581 0.114
ema_dayEMA 7. -0.200 0.145 -1.376 0.169
physical_excercise_dur 0.233 0.114 2.048 0.041
acc_delta 0.103 0.005 19.839 0.000
temp_mean_z 0.073 0.023 3.140 0.002
temp_slope_z 0.052 0.023 2.216 0.027
ema_beep_f2:ema_dayEMA 2. 0.061 0.188 0.325 0.745
ema_beep_f3:ema_dayEMA 2. 0.257 0.192 1.341 0.180
ema_beep_f4:ema_dayEMA 2. 0.110 0.189 0.584 0.559
ema_beep_f5:ema_dayEMA 2. 0.120 0.191 0.629 0.530
ema_beep_f6:ema_dayEMA 2. -0.134 0.187 -0.716 0.474
ema_beep_f2:ema_dayEMA 3. 0.078 0.191 0.410 0.682
ema_beep_f3:ema_dayEMA 3. 0.151 0.190 0.791 0.429
ema_beep_f4:ema_dayEMA 3. 0.217 0.189 1.150 0.250
ema_beep_f5:ema_dayEMA 3. 0.185 0.194 0.955 0.340
ema_beep_f6:ema_dayEMA 3. 0.249 0.192 1.297 0.195
ema_beep_f2:ema_dayEMA 4. 0.076 0.193 0.395 0.693
ema_beep_f3:ema_dayEMA 4. 0.066 0.193 0.344 0.731
ema_beep_f4:ema_dayEMA 4. -0.003 0.190 -0.014 0.989
ema_beep_f5:ema_dayEMA 4. 0.090 0.193 0.468 0.640
ema_beep_f6:ema_dayEMA 4. -0.141 0.189 -0.750 0.453
ema_beep_f2:ema_dayEMA 5. 0.017 0.187 0.090 0.928
ema_beep_f3:ema_dayEMA 5. 0.098 0.187 0.521 0.602
ema_beep_f4:ema_dayEMA 5. 0.236 0.189 1.251 0.211
ema_beep_f5:ema_dayEMA 5. 0.165 0.198 0.836 0.403
ema_beep_f6:ema_dayEMA 5. -0.088 0.187 -0.470 0.638
ema_beep_f2:ema_dayEMA 6. 0.197 0.193 1.025 0.305
ema_beep_f3:ema_dayEMA 6. 0.158 0.191 0.826 0.409
ema_beep_f4:ema_dayEMA 6. 0.140 0.192 0.730 0.465
ema_beep_f5:ema_dayEMA 6. 0.250 0.199 1.254 0.210
ema_beep_f6:ema_dayEMA 6. -0.029 0.192 -0.150 0.880
ema_beep_f2:ema_dayEMA 7. 0.276 0.198 1.392 0.164
ema_beep_f3:ema_dayEMA 7. 0.148 0.194 0.764 0.445
ema_beep_f4:ema_dayEMA 7. 0.194 0.198 0.981 0.327
ema_beep_f5:ema_dayEMA 7. 0.100 0.199 0.502 0.616
ema_beep_f6:ema_dayEMA 7. 0.029 0.214 0.136 0.892
-----------------------------------------------------------------
[4m
RANDOM EFFECTS:
[24m--------------------------------------------------------------------------
Group Parameter Std. Dev.
------------------------------------- ------------------------ -----------
castor_record_id.ema_day.ema_beep_f (Intercept) 0.527
castor_record_id acc_delta 0.032
castor_record_id.1 physical_excercise_dur 0.493
castor_record_id.2 temp_mean_z 0.134
castor_record_id.3 temp_slope_z 0.122
castor_record_id.4 (Intercept) 0.284
castor_record_id.4 Week_TypeStress 0.305
Residual 0.337
--------------------------------------------------------------------------
[4m
Grouping variables:
[24m--------------------------------------------------------
Group # groups ICC
------------------------------------- ---------- -------
castor_record_id:ema_day:ema_beep_f 2879 0.371
castor_record_id 82 0.001
--------------------------------------------------------
[34m# Check for Multicollinearity
[39m
[32mLow Correlation
[39m Parameter VIF Increased SE
Week_Type 1.00 1.00
Sex 1.02 1.01
Program 1.02 1.01
First_Week 1.00 1.00
physical_excercise_dur 1.01 1.00
acc_delta 1.02 1.01
temp_mean_z 1.04 1.02
temp_slope_z 1.01 1.01
[31mHigh Correlation
[39m Parameter VIF Increased SE
ema_beep_f 15343.87 123.87
ema_day 69541.44 263.71
ema_beep_f:ema_day 275633465.37 16602.21
# Model
# glmer.sc_ph_auc_week <- glmer( sc_phasic_auc_s ~ Week_Type +
# Sex + Program + # Model pop differences
# First_Week + ema_beep_f*ema_day + # Model day related differencxes
# physical_excercise_dur +acc_delta + # Modle movement stuff
# temp_mean_z + temp_slope_z + # Model temp effects
# (1+Week_Type |castor_record_id) +
# (0+temp_slope_z|castor_record_id) + (0+temp_mean_z|castor_record_id) +
# (0+physical_excercise_dur|castor_record_id) + (0+acc_delta|castor_record_id) +
# (1 | castor_record_id:ema_day:ema_beep_f),
# EMA_Data,
# family=Gamma(link=log),
# control=glmerControl(calc.derivs = FALSE, optCtrl=list(maxfun=100000)))
# Model Summary
model_output(glmer.sc_ph_auc_week)
Could not calculate r-squared. Try removing missing data
before fitting the model.
[4mMODEL INFO:[24m
[3mObservations:[23m 4565
[3mDependent Variable:[23m sc_phasic_auc_s
[3mType:[23m Mixed effects generalized linear regression
[3mError Distribution: [23mGamma
[3mLink function: [23mlog
[4mMODEL FIT:[24m
[3mAIC[23m = 8232.712, [3mBIC[23m = 8611.856
[4mFIXED EFFECTS:
[24m-----------------------------------------------------------------
Est. S.E. t val. p
------------------------------- -------- ------- -------- -------
(Intercept) 0.361 0.063 5.703 0.000
Week_TypeStress -0.037 0.021 -1.747 0.081
SexMale 0.024 0.034 0.713 0.476
ProgramBSc_Med 0.040 0.038 1.054 0.292
First_WeekExam-First -0.037 0.035 -1.037 0.300
ema_beep_f2 -0.081 0.070 -1.150 0.250
ema_beep_f3 -0.054 0.070 -0.774 0.439
ema_beep_f4 -0.191 0.070 -2.733 0.006
ema_beep_f5 -0.141 0.073 -1.926 0.054
ema_beep_f6 -0.036 0.071 -0.509 0.611
ema_dayEMA 2. -0.027 0.071 -0.384 0.701
ema_dayEMA 3. -0.159 0.073 -2.188 0.029
ema_dayEMA 4. -0.195 0.073 -2.685 0.007
ema_dayEMA 5. -0.197 0.072 -2.734 0.006
ema_dayEMA 6. -0.183 0.074 -2.469 0.014
ema_dayEMA 7. -0.172 0.075 -2.289 0.022
physical_excercise_dur 0.080 0.060 1.327 0.184
acc_delta 0.039 0.002 15.712 0.000
temp_mean_z 0.061 0.014 4.349 0.000
temp_slope_z 0.042 0.014 2.920 0.003
ema_beep_f2:ema_dayEMA 2. -0.017 0.098 -0.170 0.865
ema_beep_f3:ema_dayEMA 2. -0.038 0.098 -0.383 0.702
ema_beep_f4:ema_dayEMA 2. 0.070 0.098 0.715 0.475
ema_beep_f5:ema_dayEMA 2. 0.024 0.101 0.236 0.813
ema_beep_f6:ema_dayEMA 2. -0.171 0.099 -1.725 0.084
ema_beep_f2:ema_dayEMA 3. 0.100 0.100 0.995 0.320
ema_beep_f3:ema_dayEMA 3. 0.048 0.100 0.485 0.628
ema_beep_f4:ema_dayEMA 3. 0.198 0.099 2.005 0.045
ema_beep_f5:ema_dayEMA 3. 0.150 0.103 1.457 0.145
ema_beep_f6:ema_dayEMA 3. 0.133 0.101 1.323 0.186
ema_beep_f2:ema_dayEMA 4. 0.172 0.101 1.703 0.089
ema_beep_f3:ema_dayEMA 4. 0.071 0.100 0.711 0.477
ema_beep_f4:ema_dayEMA 4. 0.202 0.100 2.020 0.043
ema_beep_f5:ema_dayEMA 4. 0.144 0.102 1.413 0.158
ema_beep_f6:ema_dayEMA 4. 0.037 0.100 0.371 0.711
ema_beep_f2:ema_dayEMA 5. 0.080 0.100 0.806 0.420
ema_beep_f3:ema_dayEMA 5. 0.063 0.099 0.630 0.529
ema_beep_f4:ema_dayEMA 5. 0.247 0.099 2.481 0.013
ema_beep_f5:ema_dayEMA 5. 0.069 0.106 0.649 0.516
ema_beep_f6:ema_dayEMA 5. 0.045 0.101 0.444 0.657
ema_beep_f2:ema_dayEMA 6. 0.148 0.101 1.470 0.142
ema_beep_f3:ema_dayEMA 6. 0.074 0.100 0.738 0.461
ema_beep_f4:ema_dayEMA 6. 0.178 0.101 1.759 0.079
ema_beep_f5:ema_dayEMA 6. 0.124 0.106 1.172 0.241
ema_beep_f6:ema_dayEMA 6. 0.057 0.103 0.557 0.577
ema_beep_f2:ema_dayEMA 7. 0.130 0.103 1.261 0.207
ema_beep_f3:ema_dayEMA 7. 0.105 0.101 1.040 0.299
ema_beep_f4:ema_dayEMA 7. 0.192 0.103 1.855 0.064
ema_beep_f5:ema_dayEMA 7. 0.040 0.105 0.385 0.700
ema_beep_f6:ema_dayEMA 7. 0.046 0.114 0.401 0.688
-----------------------------------------------------------------
[4m
RANDOM EFFECTS:
[24m--------------------------------------------------------------------------
Group Parameter Std. Dev.
------------------------------------- ------------------------ -----------
castor_record_id.ema_day.ema_beep_f (Intercept) 0.305
castor_record_id acc_delta 0.017
castor_record_id.1 physical_excercise_dur 0.284
castor_record_id.2 temp_mean_z 0.094
castor_record_id.3 temp_slope_z 0.094
castor_record_id.4 (Intercept) 0.148
castor_record_id.4 Week_TypeStress 0.157
Residual 0.339
--------------------------------------------------------------------------
[4m
Grouping variables:
[24m--------------------------------------------------------
Group # groups ICC
------------------------------------- ---------- -------
castor_record_id:ema_day:ema_beep_f 2879 0.284
castor_record_id 82 0.001
--------------------------------------------------------
Model has interaction terms. VIFs might be inflated. You may check multicollinearity among predictors of a model without interaction terms.
[34m# Check for Multicollinearity
[39m
[32mLow Correlation
[39m Parameter VIF Increased SE
Week_Type 1.00 1.00
Sex 1.02 1.01
Program 1.02 1.01
First_Week 1.00 1.00
physical_excercise_dur 1.01 1.00
acc_delta 1.02 1.01
temp_mean_z 1.03 1.01
temp_slope_z 1.01 1.00
[31mHigh Correlation
[39m Parameter VIF Increased SE
ema_beep_f 11818.47 108.71
ema_day 55741.68 236.10
ema_beep_f:ema_day 181094981.80 13457.15
Next We check the heart rate data. It again looks like the relationship of the models is not as expected.
# Model
# glmer.hr_mean_week <- glmer( hr_mean_s ~ Week_Type +
# Sex + Program + # Model pop differences
# First_Week + ema_beep_f*ema_day + # Model day related differencxes
# physical_excercise_dur +acc_delta + # Modle movement stuff
# temp_mean_z + temp_slope_z + # Model temp effects
# (1+Week_Type |castor_record_id) +
# (0+temp_slope_z|castor_record_id) + (0+temp_mean_z|castor_record_id) +
# (0+physical_excercise_dur|castor_record_id) + (0+acc_delta|castor_record_id) +
# (1 | castor_record_id:ema_day:ema_beep_f),
# EMA_Data,
# family=Gamma(link=identity),
# control=glmerControl(calc.derivs = FALSE, optCtrl=list(maxfun=100000)))
# Model Summary
model_output(glmer.hr_mean_week)
Could not calculate r-squared. Try removing missing data
before fitting the model.
[4mMODEL INFO:[24m
[3mObservations:[23m 4419
[3mDependent Variable:[23m hr_mean_s
[3mType:[23m Mixed effects generalized linear regression
[3mError Distribution: [23mGamma
[3mLink function: [23midentity
[4mMODEL FIT:[24m
[3mAIC[23m = 7074.772, [3mBIC[23m = 7451.999
[4mFIXED EFFECTS:
[24m-----------------------------------------------------------------
Est. S.E. t val. p
------------------------------- -------- ------- -------- -------
(Intercept) 3.228 0.095 34.100 0.000
Week_TypeStress -0.048 0.027 -1.775 0.076
SexMale -0.045 0.058 -0.765 0.444
ProgramBSc_Med 0.152 0.064 2.390 0.017
First_WeekExam-First 0.095 0.060 1.586 0.113
ema_beep_f2 -0.003 0.098 -0.032 0.974
ema_beep_f3 -0.030 0.098 -0.302 0.763
ema_beep_f4 0.084 0.099 0.845 0.398
ema_beep_f5 0.033 0.104 0.317 0.751
ema_beep_f6 -0.326 0.096 -3.388 0.001
ema_dayEMA 2. -0.001 0.099 -0.010 0.992
ema_dayEMA 3. -0.055 0.104 -0.530 0.596
ema_dayEMA 4. 0.010 0.103 0.096 0.924
ema_dayEMA 5. -0.041 0.100 -0.407 0.684
ema_dayEMA 6. 0.035 0.105 0.333 0.739
ema_dayEMA 7. 0.006 0.106 0.054 0.957
physical_excercise_dur 0.463 0.088 5.242 0.000
acc_delta 0.126 0.004 28.724 0.000
temp_mean_z -0.057 0.023 -2.531 0.011
temp_slope_z 0.026 0.020 1.261 0.207
ema_beep_f2:ema_dayEMA 2. 0.063 0.137 0.455 0.649
ema_beep_f3:ema_dayEMA 2. 0.048 0.139 0.343 0.731
ema_beep_f4:ema_dayEMA 2. 0.072 0.140 0.511 0.609
ema_beep_f5:ema_dayEMA 2. -0.022 0.141 -0.158 0.874
ema_beep_f6:ema_dayEMA 2. 0.058 0.133 0.432 0.666
ema_beep_f2:ema_dayEMA 3. 0.073 0.142 0.510 0.610
ema_beep_f3:ema_dayEMA 3. 0.062 0.141 0.436 0.663
ema_beep_f4:ema_dayEMA 3. 0.046 0.142 0.325 0.745
ema_beep_f5:ema_dayEMA 3. 0.029 0.146 0.199 0.842
ema_beep_f6:ema_dayEMA 3. 0.041 0.137 0.301 0.763
ema_beep_f2:ema_dayEMA 4. -0.053 0.143 -0.374 0.708
ema_beep_f3:ema_dayEMA 4. 0.056 0.143 0.392 0.695
ema_beep_f4:ema_dayEMA 4. -0.090 0.143 -0.629 0.529
ema_beep_f5:ema_dayEMA 4. -0.007 0.145 -0.046 0.963
ema_beep_f6:ema_dayEMA 4. 0.042 0.136 0.307 0.759
ema_beep_f2:ema_dayEMA 5. 0.040 0.138 0.290 0.772
ema_beep_f3:ema_dayEMA 5. 0.069 0.140 0.495 0.621
ema_beep_f4:ema_dayEMA 5. -0.026 0.140 -0.184 0.854
ema_beep_f5:ema_dayEMA 5. -0.012 0.147 -0.085 0.933
ema_beep_f6:ema_dayEMA 5. -0.061 0.135 -0.455 0.649
ema_beep_f2:ema_dayEMA 6. 0.087 0.142 0.615 0.539
ema_beep_f3:ema_dayEMA 6. -0.061 0.142 -0.431 0.667
ema_beep_f4:ema_dayEMA 6. -0.026 0.144 -0.179 0.858
ema_beep_f5:ema_dayEMA 6. -0.134 0.149 -0.896 0.370
ema_beep_f6:ema_dayEMA 6. -0.115 0.138 -0.828 0.408
ema_beep_f2:ema_dayEMA 7. -0.014 0.144 -0.094 0.925
ema_beep_f3:ema_dayEMA 7. 0.016 0.143 0.112 0.911
ema_beep_f4:ema_dayEMA 7. -0.055 0.148 -0.373 0.709
ema_beep_f5:ema_dayEMA 7. 0.003 0.148 0.020 0.984
ema_beep_f6:ema_dayEMA 7. 0.175 0.155 1.128 0.259
-----------------------------------------------------------------
[4m
RANDOM EFFECTS:
[24m--------------------------------------------------------------------------
Group Parameter Std. Dev.
------------------------------------- ------------------------ -----------
castor_record_id.ema_day.ema_beep_f (Intercept) 0.402
castor_record_id acc_delta 0.029
castor_record_id.1 physical_excercise_dur 0.414
castor_record_id.2 temp_mean_z 0.158
castor_record_id.3 temp_slope_z 0.119
castor_record_id.4 (Intercept) 0.247
castor_record_id.4 Week_TypeStress 0.191
Residual 0.139
--------------------------------------------------------------------------
[4m
Grouping variables:
[24m--------------------------------------------------------
Group # groups ICC
------------------------------------- ---------- -------
castor_record_id:ema_day:ema_beep_f 2845 0.356
castor_record_id 82 0.002
--------------------------------------------------------
Model has interaction terms. VIFs might be inflated. You may check multicollinearity among predictors of a model without interaction terms.
[34m# Check for Multicollinearity
[39m
[32mLow Correlation
[39m Parameter VIF Increased SE
Week_Type 1.00 1.00
Sex 1.02 1.01
Program 1.02 1.01
First_Week 1.00 1.00
physical_excercise_dur 1.01 1.00
acc_delta 1.02 1.01
temp_mean_z 1.03 1.02
temp_slope_z 1.01 1.00
[31mHigh Correlation
[39m Parameter VIF Increased SE
ema_beep_f 12046.78 109.76
ema_day 68217.13 261.18
ema_beep_f:ema_day 224875873.84 14995.86
# Model
# glmer.hr_max_week <- glmer( hr_max_s ~ Week_Type +
# Sex + Program + # Model pop differences
# First_Week + ema_beep_f*ema_day + # Model day related differencxes
# physical_excercise_dur +acc_delta + # Modle movement stuff
# temp_mean_z + temp_slope_z + # Model temp effects
# (1+Week_Type |castor_record_id) +
# (0+temp_slope_z|castor_record_id) + (0+temp_mean_z|castor_record_id) +
# (0+physical_excercise_dur|castor_record_id) + (0+acc_delta|castor_record_id) +
# (1 | castor_record_id:ema_day:ema_beep_f),
# EMA_Data,
# family=Gamma(link="identity"),
# control=glmerControl(calc.derivs = FALSE, optCtrl=list(maxfun=100000)))
# Model SUmmary
model_output(glmer.hr_max_week)
Could not calculate r-squared. Try removing missing data
before fitting the model.
[4mMODEL INFO:[24m
[3mObservations:[23m 4419
[3mDependent Variable:[23m hr_max_s
[3mType:[23m Mixed effects generalized linear regression
[3mError Distribution: [23mGamma
[3mLink function: [23midentity
[4mMODEL FIT:[24m
[3mAIC[23m = 10351.346, [3mBIC[23m = 10728.573
[4mFIXED EFFECTS:
[24m-----------------------------------------------------------------
Est. S.E. t val. p
------------------------------- -------- ------- -------- -------
(Intercept) 3.613 0.123 29.474 0.000
Week_TypeStress -0.093 0.032 -2.909 0.004
SexMale 0.051 0.061 0.848 0.396
ProgramBSc_Med 0.102 0.066 1.551 0.121
First_WeekExam-First 0.120 0.062 1.947 0.051
ema_beep_f2 -0.072 0.142 -0.508 0.612
ema_beep_f3 -0.086 0.142 -0.608 0.543
ema_beep_f4 0.048 0.143 0.336 0.737
ema_beep_f5 -0.009 0.150 -0.058 0.954
ema_beep_f6 -0.451 0.138 -3.279 0.001
ema_dayEMA 2. -0.092 0.142 -0.643 0.520
ema_dayEMA 3. -0.099 0.150 -0.658 0.511
ema_dayEMA 4. 0.061 0.150 0.411 0.681
ema_dayEMA 5. -0.245 0.144 -1.709 0.087
ema_dayEMA 6. 0.041 0.151 0.272 0.785
ema_dayEMA 7. -0.022 0.154 -0.140 0.888
physical_excercise_dur 0.236 0.101 2.344 0.019
acc_delta 0.155 0.006 27.809 0.000
temp_mean_z -0.197 0.028 -7.062 0.000
temp_slope_z 0.018 0.027 0.659 0.510
ema_beep_f2:ema_dayEMA 2. 0.166 0.198 0.839 0.402
ema_beep_f3:ema_dayEMA 2. 0.212 0.201 1.057 0.291
ema_beep_f4:ema_dayEMA 2. 0.048 0.201 0.237 0.812
ema_beep_f5:ema_dayEMA 2. 0.087 0.203 0.428 0.669
ema_beep_f6:ema_dayEMA 2. 0.260 0.191 1.361 0.173
ema_beep_f2:ema_dayEMA 3. 0.147 0.205 0.718 0.473
ema_beep_f3:ema_dayEMA 3. 0.097 0.203 0.476 0.634
ema_beep_f4:ema_dayEMA 3. 0.113 0.205 0.551 0.582
ema_beep_f5:ema_dayEMA 3. 0.040 0.210 0.191 0.848
ema_beep_f6:ema_dayEMA 3. 0.202 0.197 1.029 0.304
ema_beep_f2:ema_dayEMA 4. -0.070 0.206 -0.341 0.733
ema_beep_f3:ema_dayEMA 4. 0.056 0.207 0.271 0.786
ema_beep_f4:ema_dayEMA 4. -0.238 0.206 -1.153 0.249
ema_beep_f5:ema_dayEMA 4. -0.102 0.209 -0.486 0.627
ema_beep_f6:ema_dayEMA 4. 0.117 0.196 0.597 0.551
ema_beep_f2:ema_dayEMA 5. 0.242 0.199 1.219 0.223
ema_beep_f3:ema_dayEMA 5. 0.293 0.201 1.460 0.144
ema_beep_f4:ema_dayEMA 5. 0.114 0.201 0.567 0.571
ema_beep_f5:ema_dayEMA 5. 0.133 0.210 0.634 0.526
ema_beep_f6:ema_dayEMA 5. 0.095 0.191 0.494 0.621
ema_beep_f2:ema_dayEMA 6. 0.025 0.205 0.124 0.901
ema_beep_f3:ema_dayEMA 6. -0.054 0.205 -0.264 0.792
ema_beep_f4:ema_dayEMA 6. -0.103 0.208 -0.494 0.622
ema_beep_f5:ema_dayEMA 6. -0.223 0.215 -1.038 0.299
ema_beep_f6:ema_dayEMA 6. 0.025 0.199 0.123 0.902
ema_beep_f2:ema_dayEMA 7. 0.057 0.209 0.272 0.785
ema_beep_f3:ema_dayEMA 7. 0.131 0.207 0.631 0.528
ema_beep_f4:ema_dayEMA 7. -0.097 0.213 -0.454 0.650
ema_beep_f5:ema_dayEMA 7. -0.070 0.213 -0.327 0.744
ema_beep_f6:ema_dayEMA 7. 0.348 0.223 1.563 0.118
-----------------------------------------------------------------
[4m
RANDOM EFFECTS:
[24m--------------------------------------------------------------------------
Group Parameter Std. Dev.
------------------------------------- ------------------------ -----------
castor_record_id.ema_day.ema_beep_f (Intercept) 0.578
castor_record_id acc_delta 0.033
castor_record_id.1 physical_excercise_dur 0.222
castor_record_id.2 temp_mean_z 0.179
castor_record_id.3 temp_slope_z 0.137
castor_record_id.4 (Intercept) 0.236
castor_record_id.4 Week_TypeStress 0.193
Residual 0.177
--------------------------------------------------------------------------
[4m
Grouping variables:
[24m--------------------------------------------------------
Group # groups ICC
------------------------------------- ---------- -------
castor_record_id:ema_day:ema_beep_f 2845 0.639
castor_record_id 82 0.002
--------------------------------------------------------
Model has interaction terms. VIFs might be inflated. You may check multicollinearity among predictors of a model without interaction terms.
[34m# Check for Multicollinearity
[39m
[32mLow Correlation
[39m Parameter VIF Increased SE
Week_Type 1.01 1.00
Sex 1.02 1.01
Program 1.02 1.01
First_Week 1.01 1.00
physical_excercise_dur 1.01 1.01
acc_delta 1.03 1.01
temp_mean_z 1.04 1.02
temp_slope_z 1.01 1.01
[31mHigh Correlation
[39m Parameter VIF Increased SE
ema_beep_f 12191.54 110.42
ema_day 73391.98 270.91
ema_beep_f:ema_day 241939982.63 15554.42
# Model
# glmer.hr_min_week <- glmer( hr_min_s ~ Week_Type +
# Sex + Program + # Model pop differences
# First_Week + ema_beep_f*ema_day + # Model day related differencxes
# physical_excercise_dur +acc_delta + # Modle movement stuff
# temp_mean_z + temp_slope_z + # Model temp effects
# (1+Week_Type |castor_record_id) +
# (0+temp_slope_z|castor_record_id) + (0+temp_mean_z|castor_record_id) +
# (0+physical_excercise_dur|castor_record_id) + (0+acc_delta|castor_record_id) +
# (1 | castor_record_id:ema_day:ema_beep_f),
# EMA_Data,
# family=Gamma(link="identity"),
# control=glmerControl(calc.derivs = FALSE, optCtrl=list(maxfun=100000)))
# Model Summary
model_output(glmer.hr_min_week)
Could not calculate r-squared. Try removing missing data
before fitting the model.
[4mMODEL INFO:[24m
[3mObservations:[23m 4419
[3mDependent Variable:[23m hr_min_s
[3mType:[23m Mixed effects generalized linear regression
[3mError Distribution: [23mGamma
[3mLink function: [23midentity
[4mMODEL FIT:[24m
[3mAIC[23m = 9024.946, [3mBIC[23m = 9402.172
[4mFIXED EFFECTS:
[24m-----------------------------------------------------------------
Est. S.E. t val. p
------------------------------- -------- ------- -------- -------
(Intercept) 3.803 0.117 32.502 0.000
Week_TypeStress -0.034 0.035 -0.973 0.331
SexMale -0.333 0.072 -4.658 0.000
ProgramBSc_Med 0.027 0.078 0.343 0.731
First_WeekExam-First 0.039 0.073 0.527 0.598
ema_beep_f2 -0.027 0.123 -0.216 0.829
ema_beep_f3 -0.042 0.123 -0.343 0.732
ema_beep_f4 0.140 0.124 1.134 0.257
ema_beep_f5 0.013 0.130 0.098 0.922
ema_beep_f6 -0.261 0.121 -2.152 0.031
ema_dayEMA 2. -0.045 0.123 -0.364 0.716
ema_dayEMA 3. -0.047 0.129 -0.363 0.717
ema_dayEMA 4. -0.094 0.128 -0.738 0.461
ema_dayEMA 5. 0.011 0.125 0.089 0.929
ema_dayEMA 6. 0.085 0.131 0.647 0.517
ema_dayEMA 7. -0.094 0.131 -0.715 0.475
physical_excercise_dur 0.527 0.106 4.987 0.000
acc_delta 0.088 0.004 20.052 0.000
temp_mean_z 0.073 0.025 2.876 0.004
temp_slope_z 0.038 0.025 1.511 0.131
ema_beep_f2:ema_dayEMA 2. 0.099 0.172 0.575 0.565
ema_beep_f3:ema_dayEMA 2. 0.029 0.173 0.168 0.867
ema_beep_f4:ema_dayEMA 2. 0.083 0.175 0.475 0.635
ema_beep_f5:ema_dayEMA 2. 0.009 0.176 0.054 0.957
ema_beep_f6:ema_dayEMA 2. 0.053 0.168 0.313 0.754
ema_beep_f2:ema_dayEMA 3. 0.051 0.178 0.288 0.774
ema_beep_f3:ema_dayEMA 3. 0.035 0.176 0.200 0.841
ema_beep_f4:ema_dayEMA 3. -0.045 0.177 -0.256 0.798
ema_beep_f5:ema_dayEMA 3. 0.084 0.183 0.460 0.646
ema_beep_f6:ema_dayEMA 3. -0.020 0.173 -0.117 0.907
ema_beep_f2:ema_dayEMA 4. 0.014 0.177 0.077 0.939
ema_beep_f3:ema_dayEMA 4. 0.046 0.177 0.257 0.797
ema_beep_f4:ema_dayEMA 4. -0.022 0.178 -0.122 0.903
ema_beep_f5:ema_dayEMA 4. 0.205 0.181 1.131 0.258
ema_beep_f6:ema_dayEMA 4. 0.116 0.172 0.679 0.497
ema_beep_f2:ema_dayEMA 5. -0.016 0.173 -0.094 0.925
ema_beep_f3:ema_dayEMA 5. -0.038 0.175 -0.217 0.828
ema_beep_f4:ema_dayEMA 5. -0.156 0.175 -0.890 0.373
ema_beep_f5:ema_dayEMA 5. -0.034 0.184 -0.188 0.851
ema_beep_f6:ema_dayEMA 5. -0.089 0.171 -0.519 0.604
ema_beep_f2:ema_dayEMA 6. 0.029 0.178 0.163 0.870
ema_beep_f3:ema_dayEMA 6. -0.074 0.177 -0.417 0.677
ema_beep_f4:ema_dayEMA 6. -0.181 0.180 -1.004 0.315
ema_beep_f5:ema_dayEMA 6. -0.151 0.187 -0.807 0.420
ema_beep_f6:ema_dayEMA 6. -0.208 0.175 -1.188 0.235
ema_beep_f2:ema_dayEMA 7. 0.087 0.179 0.488 0.626
ema_beep_f3:ema_dayEMA 7. 0.066 0.177 0.374 0.708
ema_beep_f4:ema_dayEMA 7. -0.021 0.183 -0.113 0.910
ema_beep_f5:ema_dayEMA 7. 0.050 0.184 0.272 0.785
ema_beep_f6:ema_dayEMA 7. 0.299 0.195 1.534 0.125
-----------------------------------------------------------------
[4m
RANDOM EFFECTS:
[24m--------------------------------------------------------------------------
Group Parameter Std. Dev.
------------------------------------- ------------------------ -----------
castor_record_id.ema_day.ema_beep_f (Intercept) 0.508
castor_record_id acc_delta 0.026
castor_record_id.1 physical_excercise_dur 0.461
castor_record_id.2 temp_mean_z 0.165
castor_record_id.3 temp_slope_z 0.151
castor_record_id.4 (Intercept) 0.296
castor_record_id.4 Week_TypeStress 0.256
Residual 0.167
--------------------------------------------------------------------------
[4m
Grouping variables:
[24m--------------------------------------------------------
Group # groups ICC
------------------------------------- ---------- -------
castor_record_id:ema_day:ema_beep_f 2845 0.406
castor_record_id 82 0.001
--------------------------------------------------------
Model has interaction terms. VIFs might be inflated. You may check multicollinearity among predictors of a model without interaction terms.
[34m# Check for Multicollinearity
[39m
[32mLow Correlation
[39m Parameter VIF Increased SE
Week_Type 1.00 1.00
Sex 1.02 1.01
Program 1.02 1.01
First_Week 1.00 1.00
physical_excercise_dur 1.01 1.00
acc_delta 1.02 1.01
temp_mean_z 1.04 1.02
temp_slope_z 1.01 1.00
[31mHigh Correlation
[39m Parameter VIF Increased SE
ema_beep_f 12502.42 111.81
ema_day 60554.17 246.08
ema_beep_f:ema_day 207552984.59 14406.70
We now present the physiology results in a table, corrected for multiple comparisons with FDR.
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
htmlTable::htmlTable(tab_phy_v_week$knitr)
There were 27 warnings (use warnings() to see them)
Finally, we can plot the estimates for the week differences here, also correcting for multiple comparisons with FDR.
# Together
plot.physio <- plot_models(glmer.posmood_week, glmer.sc_ton_mean_week, glmer.sc_ph_num)
rm_term <- as.vector(plot.physio$data$term[1:97])
rm_term <- rm_term[rm_term != "Week_TypeStress"]
# Color list
col_list <-c('#7570B3','#7570B3','#7570B3',"#D95F02", "#D95F02", '#D95F02', '#D95F02', "#1B9E77", "#1B9E77")
# Make the plot
plot.all_week_resid <- plot_models( glmer.posmood_week, glmer.negmood_week,
glmer.sc_ton_mean_week, glmer.sc_ph_num, glmer.sc_ph_mag, glmer.sc_ph_auc_week,
glmer.hr_mean_week, glmer.hr_min_week, glmer.hr_max_week,
m.labels=c("Positive Affect", "Negative Affect", "SC Tonic", "SC Number", "SC Magnitued", "SC AUC","HR Mean", "HR Minimum", "HR Maximum"),
axis.labels = c(" "),
# Stastistical Stuff
rm.terms = rm_term,
#show.values = T,
#value.size = 4,
#std.est=T,
show.p=T,
p.shape=T,
p.adjust = "fdr",
legend.pval.title = "Significance",
# Visual Stuff
colors=col_list,
dot.size=2,
line.size = 1,
spacing=0.7,
vline.color = "darkgrey",
legend.title = "") +
ptheme +
theme(axis.ticks.y=element_blank(),
legend.text=element_text(size=16), legend.title = element_text(size=16),
axis.title = element_text(size=16),axis.text.x=element_text(size=11),
panel.grid.major.x = element_line(colour = "grey95"),
panel.grid.minor.x = element_line(colour = "grey90"))
# Plot and save
plot.all_week_resid + coord_flip(ylim=c(-0.5,0.5)) + theme(panel.grid.major.y = element_line(colour = "grey95"), panel.grid.minor.y = element_line(colour = "grey90"))
Coordinate system already present. Adding new coordinate system, which will replace the existing one.
ggsave("figures/fig_WeekResid_ALL.tiff", units="in", width=4, height=4, dpi=300, compression = 'lzw', bg="transparent")
TIFFOpen: figures/fig_WeekResid_ALL.tiff: No such file or directory.
We saw lower physiological arousal in the exam week, which was unexpected. We also see lower positive arousal. So it could be that the effects we see are actually driven by our positive arousal/mood. To check this, and to also check whether we can replicate previous indings on moment-to-moment associations between subjective stress and physiological arousal and mood, we discard week type, and look at the measure continuously.
In this models, we use a fixed slope due to convergence issues. We initially tried changing optimizers, which did not result in better convergence, so we simplified the model by doing a fixed slope for our fixed effect of interest. We again use more or less the same covariates as before, except for program, and for day. We keep the time variable in there to control for circadian rhythms. Additionally, we see there is a correlation between our subjective stress measures, so we need to factor this effect in by modeling interaction terms as well.
rcorr(as.matrix(EMA_Data[,c("event_tot_s", "activity_tot_s", "social_tot_s", "physical_tot_s")]))
event_tot_s activity_tot_s social_tot_s physical_tot_s
event_tot_s 1.00 0.42 0.17 0.07
activity_tot_s 0.42 1.00 0.34 0.16
social_tot_s 0.17 0.34 1.00 0.18
physical_tot_s 0.07 0.16 0.18 1.00
n= 5735
P
event_tot_s activity_tot_s social_tot_s physical_tot_s
event_tot_s 0 0 0
activity_tot_s 0 0 0
social_tot_s 0 0 0
physical_tot_s 0 0 0
First we check the moment-to-moment associations between positive affect and subjective stress. These confirm that momentary subjective stress is also associated with reduced positive affect, and increased negative affect.
# glmer.posmood_beep <- lmer(mood_positive_s ~
# event_tot_s*activity_tot_s + activity_tot_s*social_tot_s + physical_tot_s +
# ema_beep + # I will model the beep to factor circadian rhythms
# Sex + physical_excercise_dur + acc_delta +
# temp_mean_z + temp_slope_z +
# (0+event_tot_s+activity_tot_s+social_tot_s+physical_tot_s | castor_record_id) +
# (0+ema_beep|castor_record_id)+
# (0+temp_slope_z|castor_record_id) + (0+temp_mean_z|castor_record_id)+
# (0+acc_delta|castor_record_id) + (0+physical_excercise_dur|castor_record_id),
# EMA_Data,
# # family=gaussian(link="identity"),
# control=lmerControl(calc.derivs = FALSE) )
# Summary
model_output(glmer.posmood_beep )
[4mMODEL INFO:[24m
[3mObservations:[23m 4793
[3mDependent Variable:[23m mood_positive_s
[3mType:[23m Mixed effects linear regression
[4mMODEL FIT:[24m
[3mAIC[23m = 14811.445, [3mBIC[23m = 14999.218
[3mPseudo-R² (fixed effects)[23m = 0.257
[3mPseudo-R² (total)[23m = 0.530
[4mFIXED EFFECTS:
[24m------------------------------------------------------------------------------
Est. S.E. t val. d.f. p
--------------------------------- -------- ------- -------- ---------- -------
(Intercept) 9.714 0.204 47.520 4193.440 0.000
event_tot_s -0.120 0.035 -3.452 664.640 0.001
activity_tot_s -0.200 0.040 -4.951 2124.591 0.000
social_tot_s -0.299 0.031 -9.745 629.127 0.000
physical_tot_s -0.152 0.020 -7.628 92.148 0.000
ema_beep -0.010 0.013 -0.816 122.407 0.416
SexMale 0.271 0.142 1.915 246.186 0.057
physical_excercise_dur 0.530 0.120 4.426 36.837 0.000
acc_delta 0.004 0.003 1.298 4597.881 0.194
temp_mean_z -0.003 0.028 -0.109 59.550 0.913
temp_slope_z -0.015 0.016 -0.936 4581.209 0.350
event_tot_s:activity_tot_s -0.017 0.006 -2.808 3581.305 0.005
activity_tot_s:social_tot_s 0.027 0.005 5.131 2882.146 0.000
------------------------------------------------------------------------------
[3m[36mp values calculated using Satterthwaite d.f.[39m[23m
[4m
RANDOM EFFECTS:
[24m---------------------------------------------------------
Group Parameter Std. Dev.
-------------------- ------------------------ -----------
castor_record_id event_tot_s 0.131
castor_record_id activity_tot_s 0.071
castor_record_id social_tot_s 0.111
castor_record_id physical_tot_s 0.135
castor_record_id.1 ema_beep 0.068
castor_record_id.2 temp_slope_z 0.000
castor_record_id.3 temp_mean_z 0.151
castor_record_id.4 acc_delta 0.000
castor_record_id.5 physical_excercise_dur 0.308
Residual 1.056
---------------------------------------------------------
[4m
Grouping variables:
[24m-------------------------------------
Group # groups ICC
------------------ ---------- -------
castor_record_id 82 0.014
-------------------------------------
Model has interaction terms. VIFs might be inflated. You may check multicollinearity among predictors of a model without interaction terms.
[34m# Check for Multicollinearity
[39m
[32mLow Correlation
[39m Parameter VIF Increased SE
event_tot_s 3.41 1.85
social_tot_s 3.40 1.84
physical_tot_s 1.19 1.09
ema_beep 1.04 1.02
Sex 1.00 1.00
physical_excercise_dur 1.01 1.00
acc_delta 1.03 1.02
temp_mean_z 1.03 1.02
temp_slope_z 1.02 1.01
[33mModerate Correlation
[39m Parameter VIF Increased SE
activity_tot_s:social_tot_s 7.92 2.81
[31mHigh Correlation
[39m Parameter VIF Increased SE
activity_tot_s 12.52 3.54
event_tot_s:activity_tot_s 16.21 4.03
# glmer.negmood_beep <- glmer(mood_negative_s ~
# event_tot_s*activity_tot_s + activity_tot_s*social_tot_s + physical_tot_s +
# ema_beep + # I will model the beep to factor circadian rhythms
# Sex + physical_excercise_dur + acc_delta +
# temp_mean_z + temp_slope_z +
# (0+event_tot_s+activity_tot_s+social_tot_s+physical_tot_s | castor_record_id) +
# (0+ema_beep|castor_record_id)+
# (0+temp_slope_z|castor_record_id) + (0+temp_mean_z|castor_record_id)+
# (0+acc_delta|castor_record_id) + (0+physical_excercise_dur|castor_record_id),
# EMA_Data,
# family=Gamma(link="identity"),
# control=glmerControl(calc.derivs = FALSE) )
# Summary
model_output(glmer.negmood_beep)
[4mMODEL INFO:[24m
[3mObservations:[23m 4793
[3mDependent Variable:[23m mood_negative_s
[3mType:[23m Mixed effects generalized linear regression
[3mError Distribution: [23mGamma
[3mLink function: [23midentity
[4mMODEL FIT:[24m
[3mAIC[23m = 11818.666, [3mBIC[23m = 12006.439
[3mPseudo-R² (fixed effects)[23m = NA
[3mPseudo-R² (total)[23m = NA
[4mFIXED EFFECTS:
[24m-------------------------------------------------------------------
Est. S.E. t val. p
--------------------------------- -------- ------- -------- -------
(Intercept) 0.775 0.134 5.788 0.000
event_tot_s 0.035 0.021 1.648 0.099
activity_tot_s 0.080 0.029 2.749 0.006
social_tot_s 0.234 0.023 10.148 0.000
physical_tot_s 0.127 0.015 8.211 0.000
ema_beep 0.020 0.010 2.024 0.043
SexMale -0.244 0.096 -2.531 0.011
physical_excercise_dur -0.195 0.080 -2.444 0.015
acc_delta -0.002 0.003 -0.705 0.481
temp_mean_z 0.021 0.025 0.823 0.410
temp_slope_z -0.013 0.015 -0.832 0.406
event_tot_s:activity_tot_s 0.011 0.004 2.612 0.009
activity_tot_s:social_tot_s -0.011 0.004 -2.712 0.007
-------------------------------------------------------------------
[4m
RANDOM EFFECTS:
[24m---------------------------------------------------------
Group Parameter Std. Dev.
-------------------- ------------------------ -----------
castor_record_id event_tot_s 0.061
castor_record_id activity_tot_s 0.066
castor_record_id social_tot_s 0.096
castor_record_id physical_tot_s 0.101
castor_record_id.1 ema_beep 0.057
castor_record_id.2 temp_slope_z 0.057
castor_record_id.3 temp_mean_z 0.155
castor_record_id.4 acc_delta 0.012
castor_record_id.5 physical_excercise_dur 0.226
Residual 0.376
---------------------------------------------------------
[4m
Grouping variables:
[24m-------------------------------------
Group # groups ICC
------------------ ---------- -------
castor_record_id 82 0.016
-------------------------------------
Model has interaction terms. VIFs might be inflated. You may check multicollinearity among predictors of a model without interaction terms.
[34m# Check for Multicollinearity
[39m
[32mLow Correlation
[39m Parameter VIF Increased SE
event_tot_s 3.21 1.79
social_tot_s 2.67 1.64
physical_tot_s 1.09 1.05
ema_beep 1.02 1.01
Sex 1.01 1.00
physical_excercise_dur 1.01 1.00
acc_delta 1.02 1.01
temp_mean_z 1.02 1.01
temp_slope_z 1.01 1.01
[33mModerate Correlation
[39m Parameter VIF Increased SE
activity_tot_s 8.16 2.86
activity_tot_s:social_tot_s 5.50 2.35
[31mHigh Correlation
[39m Parameter VIF Increased SE
event_tot_s:activity_tot_s 11.22 3.35
tab_model(glmer.posmood_beep, glmer.negmood_beep,
# terms=c("(Intercept)","Week_Type [Stress]", "Sex [Male]", "Program [BSc_Med]", "First_Week [Exam-First]", "acc_delta", "physical_excercise_dur"),
title="Table 3. HR vs Week Type",
transform=NULL, p.adjust = "fdr",
show.stat=TRUE,
show.se=TRUE) %>%
return() %$%
knitr %>%
asis_output()
| mood positive s | mood negative s | |||||||||
|---|---|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | Statistic | p | Estimates | std. Error | CI | Statistic | p |
| (Intercept) | 9.71 | 0.20 | 9.31 – 10.11 | 47.52 | <0.001 | 0.77 | 0.13 | 0.51 – 1.04 | 5.79 | <0.001 |
| event_tot_s | -0.12 | 0.03 | -0.19 – -0.05 | -3.45 | 0.001 | 0.04 | 0.02 | -0.01 – 0.08 | 1.65 | 0.129 |
| activity_tot_s | -0.20 | 0.04 | -0.28 – -0.12 | -4.95 | <0.001 | 0.08 | 0.03 | 0.02 – 0.14 | 2.75 | 0.017 |
| social_tot_s | -0.30 | 0.03 | -0.36 – -0.24 | -9.75 | <0.001 | 0.23 | 0.02 | 0.19 – 0.28 | 10.15 | <0.001 |
| physical_tot_s | -0.15 | 0.02 | -0.19 – -0.11 | -7.63 | <0.001 | 0.13 | 0.02 | 0.10 – 0.16 | 8.21 | <0.001 |
| ema_beep | -0.01 | 0.01 | -0.03 – 0.01 | -0.82 | 0.449 | 0.02 | 0.01 | 0.00 – 0.04 | 2.02 | 0.062 |
| Sex [Male] | 0.27 | 0.14 | -0.01 – 0.55 | 1.92 | 0.080 | -0.24 | 0.10 | -0.43 – -0.05 | -2.53 | 0.021 |
| physical_excercise_dur | 0.53 | 0.12 | 0.30 – 0.76 | 4.43 | <0.001 | -0.19 | 0.08 | -0.35 – -0.04 | -2.44 | 0.024 |
| acc_delta | 0.00 | 0.00 | -0.00 – 0.01 | 1.30 | 0.253 | -0.00 | 0.00 | -0.01 – 0.00 | -0.71 | 0.481 |
| temp_mean_z | -0.00 | 0.03 | -0.06 – 0.05 | -0.11 | 0.913 | 0.02 | 0.03 | -0.03 – 0.07 | 0.82 | 0.445 |
| temp_slope_z | -0.02 | 0.02 | -0.05 – 0.02 | -0.94 | 0.413 | -0.01 | 0.02 | -0.04 – 0.02 | -0.83 | 0.445 |
|
event_tot_s * activity_tot_s |
-0.02 | 0.01 | -0.03 – -0.01 | -2.81 | 0.008 | 0.01 | 0.00 | 0.00 – 0.02 | 2.61 | 0.019 |
|
activity_tot_s * social_tot_s |
0.03 | 0.01 | 0.02 – 0.04 | 5.13 | <0.001 | -0.01 | 0.00 | -0.02 – -0.00 | -2.71 | 0.017 |
| Random Effects | ||||||||||
| σ2 | 1.12 | 0.14 | ||||||||
| τ00 | 0.02 castor_record_id | 0.00 castor_record_id | ||||||||
| 0.00 castor_record_id.1 | 0.00 castor_record_id.1 | |||||||||
| 0.00 castor_record_id.2 | 0.00 castor_record_id.2 | |||||||||
| 0.02 castor_record_id.3 | 0.02 castor_record_id.3 | |||||||||
| 0.00 castor_record_id.4 | 0.00 castor_record_id.4 | |||||||||
| 0.10 castor_record_id.5 | 0.05 castor_record_id.5 | |||||||||
| τ11 | 0.01 castor_record_id.activity_tot_s | 0.00 castor_record_id.activity_tot_s | ||||||||
| 0.01 castor_record_id.social_tot_s | 0.01 castor_record_id.social_tot_s | |||||||||
| 0.02 castor_record_id.physical_tot_s | 0.01 castor_record_id.physical_tot_s | |||||||||
| ρ01 | -0.23 castor_record_id.activity_tot_s | -0.26 castor_record_id.activity_tot_s | ||||||||
| -0.27 castor_record_id.social_tot_s | -0.19 castor_record_id.social_tot_s | |||||||||
| -0.21 castor_record_id.physical_tot_s | -0.20 castor_record_id.physical_tot_s | |||||||||
| ICC | 0.74 | |||||||||
| N | 82 castor_record_id | 82 castor_record_id | ||||||||
| Observations | 4793 | 4793 | ||||||||
| Marginal R2 / Conditional R2 | 0.353 / NA | 0.357 / 0.833 | ||||||||
We next model the skin conductance vs. the subjective stress measures. We expect that continuous analysis will result in a replication of previous findings where momentary stress is indeed associated with increase physiological arousal. Our findings are indeed shown in Table: Subjective Stress and Outcomes.
# Model
# glmer.sc_ton_ema <- glmer(sc_tonic_mean_s ~
# event_tot_s*activity_tot_s + activity_tot_s*social_tot_s + physical_tot_s +
# ema_beep + # I will model the beep to factor circadian rhythms
# Sex + physical_excercise_dur + acc_delta +
# temp_mean_z + temp_slope_z +
# (0+event_tot_s+activity_tot_s+social_tot_s+physical_tot_s | castor_record_id) +
# (0+ema_beep|castor_record_id)+
# (0+temp_slope_z|castor_record_id) + (0+temp_mean_z|castor_record_id)+
# (0+acc_delta|castor_record_id) + (0+physical_excercise_dur|castor_record_id),
# EMA_Data,
# family=Gamma(link="log"),
# control=glmerControl(calc.derivs = FALSE) )
# Summary
model_output(glmer.sc_ton_ema)
[4mMODEL INFO:[24m
[3mObservations:[23m 4565
[3mDependent Variable:[23m sc_tonic_mean_s
[3mType:[23m Mixed effects generalized linear regression
[3mError Distribution: [23mGamma
[3mLink function: [23mlog
[4mMODEL FIT:[24m
[3mAIC[23m = 703.846, [3mBIC[23m = 890.205
[3mPseudo-R² (fixed effects)[23m = NA
[3mPseudo-R² (total)[23m = NA
[4mFIXED EFFECTS:
[24m-------------------------------------------------------------------
Est. S.E. t val. p
--------------------------------- -------- ------- -------- -------
(Intercept) -0.005 0.054 -0.089 0.929
event_tot_s 0.013 0.009 1.502 0.133
activity_tot_s 0.014 0.011 1.291 0.197
social_tot_s 0.000 0.008 0.009 0.993
physical_tot_s -0.002 0.004 -0.533 0.594
ema_beep 0.001 0.003 0.428 0.669
SexMale 0.001 0.020 0.060 0.952
physical_excercise_dur 0.083 0.049 1.688 0.091
acc_delta 0.018 0.002 9.518 0.000
temp_mean_z 0.042 0.009 4.487 0.000
temp_slope_z 0.027 0.009 2.835 0.005
event_tot_s:activity_tot_s -0.003 0.002 -1.643 0.100
activity_tot_s:social_tot_s -0.001 0.001 -0.445 0.656
-------------------------------------------------------------------
[4m
RANDOM EFFECTS:
[24m---------------------------------------------------------
Group Parameter Std. Dev.
-------------------- ------------------------ -----------
castor_record_id event_tot_s 0.023
castor_record_id activity_tot_s 0.020
castor_record_id social_tot_s 0.024
castor_record_id physical_tot_s 0.024
castor_record_id.1 ema_beep 0.015
castor_record_id.2 temp_slope_z 0.062
castor_record_id.3 temp_mean_z 0.065
castor_record_id.4 acc_delta 0.014
castor_record_id.5 physical_excercise_dur 0.289
Residual 0.279
---------------------------------------------------------
[4m
Grouping variables:
[24m-------------------------------------
Group # groups ICC
------------------ ---------- -------
castor_record_id 82 0.003
-------------------------------------
Model has interaction terms. VIFs might be inflated. You may check multicollinearity among predictors of a model without interaction terms.
[34m# Check for Multicollinearity
[39m
[32mLow Correlation
[39m Parameter VIF Increased SE
social_tot_s 4.73 2.17
physical_tot_s 1.11 1.05
ema_beep 1.04 1.02
Sex 1.01 1.00
physical_excercise_dur 1.00 1.00
acc_delta 1.01 1.00
temp_mean_z 1.02 1.01
temp_slope_z 1.00 1.00
[33mModerate Correlation
[39m Parameter VIF Increased SE
event_tot_s 5.32 2.31
[31mHigh Correlation
[39m Parameter VIF Increased SE
activity_tot_s 13.47 3.67
event_tot_s:activity_tot_s 21.44 4.63
activity_tot_s:social_tot_s 11.33 3.37
# Model
# glmer.sc_ph_num_ema <- glmer(sc_phasic_num ~
# event_tot_s*activity_tot_s + activity_tot_s*social_tot_s + physical_tot_s +
# ema_beep + # I will model the beep to factor circadian rhythms
# Sex + physical_excercise_dur + acc_delta +
# temp_mean_z + temp_slope_z +
# (0+event_tot_s+activity_tot_s+social_tot_s+physical_tot_s | castor_record_id) +
# (0+ema_beep|castor_record_id)+
# (0+temp_slope_z|castor_record_id) + (0+temp_mean_z|castor_record_id)+
# (0+acc_delta|castor_record_id) + (0+physical_excercise_dur|castor_record_id),
# EMA_Data,
# family=poisson(link="sqrt"),
# control=glmerControl(calc.derivs = FALSE) )
# Model Summary
model_output(glmer.sc_ph_num_ema)
Could not calculate r-squared. Try removing missing data
before fitting the model.
[4mMODEL INFO:[24m
[3mObservations:[23m 4565
[3mDependent Variable:[23m sc_phasic_num
[3mType:[23m Mixed effects generalized linear regression
[3mError Distribution: [23mpoisson
[3mLink function: [23msqrt
[4mMODEL FIT:[24m
[3mAIC[23m = 75466.975, [3mBIC[23m = 75646.908
[4mFIXED EFFECTS:
[24m-------------------------------------------------------------------
Est. S.E. z val. p
--------------------------------- -------- ------- -------- -------
(Intercept) 2.006 0.113 17.691 0.000
event_tot_s 0.039 0.029 1.349 0.177
activity_tot_s 0.083 0.038 2.201 0.028
social_tot_s 0.027 0.030 0.904 0.366
physical_tot_s 0.035 0.030 1.156 0.248
ema_beep -0.052 0.024 -2.140 0.032
SexMale 0.843 0.115 7.362 0.000
physical_excercise_dur -0.197 0.409 -0.482 0.630
acc_delta 0.228 0.014 16.039 0.000
temp_mean_z 0.427 0.091 4.710 0.000
temp_slope_z 0.223 0.069 3.236 0.001
event_tot_s:activity_tot_s -0.009 0.003 -2.782 0.005
activity_tot_s:social_tot_s -0.011 0.003 -3.916 0.000
-------------------------------------------------------------------
[4m
RANDOM EFFECTS:
[24m---------------------------------------------------------
Group Parameter Std. Dev.
-------------------- ------------------------ -----------
castor_record_id event_tot_s 0.211
castor_record_id activity_tot_s 0.279
castor_record_id social_tot_s 0.234
castor_record_id physical_tot_s 0.266
castor_record_id.1 ema_beep 0.214
castor_record_id.2 temp_slope_z 0.610
castor_record_id.3 temp_mean_z 0.810
castor_record_id.4 acc_delta 0.126
castor_record_id.5 physical_excercise_dur 3.508
---------------------------------------------------------
[4m
Grouping variables:
[24m-------------------------------------
Group # groups ICC
------------------ ---------- -------
castor_record_id 82 0.003
-------------------------------------
[34m# Check for Multicollinearity
[39m
[32mLow Correlation
[39m Parameter VIF Increased SE
event_tot_s 1.83 1.35
activity_tot_s 2.80 1.67
social_tot_s 2.22 1.49
physical_tot_s 1.31 1.14
ema_beep 1.00 1.00
Sex 1.00 1.00
physical_excercise_dur 1.00 1.00
acc_delta 1.00 1.00
temp_mean_z 1.00 1.00
temp_slope_z 1.00 1.00
event_tot_s:activity_tot_s 2.77 1.66
activity_tot_s:social_tot_s 2.32 1.52
# Model
# glmer.sc_ph_mag_ema <- glmer(sc_phasic_mag_s ~
# event_tot_s*activity_tot_s + activity_tot_s*social_tot_s + physical_tot_s +
# ema_beep + # I will model the beep to factor circadian rhythms
# Sex + physical_excercise_dur + acc_delta +
# temp_mean_z + temp_slope_z +
# (0+event_tot_s+activity_tot_s+social_tot_s+physical_tot_s | castor_record_id) +
# (0+ema_beep|castor_record_id)+
# (0+temp_slope_z|castor_record_id) + (0+temp_mean_z|castor_record_id)+
# (0+acc_delta|castor_record_id) + (0+physical_excercise_dur|castor_record_id),
# EMA_Data,
# family=Gamma(link="identity"),
# control=glmerControl(calc.derivs = FALSE, optCtrl=list(maxfun=100000)))
# Model Summary
model_output(glmer.sc_ph_mag_ema)
[4mMODEL INFO:[24m
[3mObservations:[23m 4565
[3mDependent Variable:[23m sc_phasic_mag_s
[3mType:[23m Mixed effects generalized linear regression
[3mError Distribution: [23mGamma
[3mLink function: [23midentity
[4mMODEL FIT:[24m
[3mAIC[23m = 11198.234, [3mBIC[23m = 11384.593
[3mPseudo-R² (fixed effects)[23m = NA
[3mPseudo-R² (total)[23m = NA
[4mFIXED EFFECTS:
[24m-------------------------------------------------------------------
Est. S.E. t val. p
--------------------------------- -------- ------- -------- -------
(Intercept) 1.338 0.165 8.124 0.000
event_tot_s 0.071 0.027 2.670 0.008
activity_tot_s 0.070 0.033 2.155 0.031
social_tot_s -0.004 0.023 -0.177 0.859
physical_tot_s 0.028 0.013 2.142 0.032
ema_beep -0.001 0.009 -0.071 0.944
SexMale 0.042 0.072 0.583 0.560
physical_excercise_dur 0.362 0.130 2.788 0.005
acc_delta 0.107 0.006 18.206 0.000
temp_mean_z 0.119 0.025 4.814 0.000
temp_slope_z 0.068 0.021 3.231 0.001
event_tot_s:activity_tot_s -0.012 0.005 -2.447 0.014
activity_tot_s:social_tot_s -0.001 0.004 -0.322 0.747
-------------------------------------------------------------------
[4m
RANDOM EFFECTS:
[24m---------------------------------------------------------
Group Parameter Std. Dev.
-------------------- ------------------------ -----------
castor_record_id event_tot_s 0.061
castor_record_id activity_tot_s 0.053
castor_record_id social_tot_s 0.050
castor_record_id physical_tot_s 0.066
castor_record_id.1 ema_beep 0.044
castor_record_id.2 temp_slope_z 0.088
castor_record_id.3 temp_mean_z 0.151
castor_record_id.4 acc_delta 0.037
castor_record_id.5 physical_excercise_dur 0.628
Residual 0.408
---------------------------------------------------------
[4m
Grouping variables:
[24m-------------------------------------
Group # groups ICC
------------------ ---------- -------
castor_record_id 82 0.006
-------------------------------------
Model has interaction terms. VIFs might be inflated. You may check multicollinearity among predictors of a model without interaction terms.
[34m# Check for Multicollinearity
[39m
[32mLow Correlation
[39m Parameter VIF Increased SE
event_tot_s 4.83 2.20
physical_tot_s 1.10 1.05
ema_beep 1.05 1.03
Sex 1.01 1.00
physical_excercise_dur 1.00 1.00
acc_delta 1.01 1.01
temp_mean_z 1.02 1.01
temp_slope_z 1.01 1.00
[33mModerate Correlation
[39m Parameter VIF Increased SE
social_tot_s 5.32 2.31
[31mHigh Correlation
[39m Parameter VIF Increased SE
activity_tot_s 14.63 3.83
event_tot_s:activity_tot_s 21.24 4.61
activity_tot_s:social_tot_s 12.72 3.57
# Model
# glmer.sc_ph_auc_ema <- glmer(sc_phasic_auc_s ~
# event_tot_s*activity_tot_s + activity_tot_s*social_tot_s + physical_tot_s +
# ema_beep + # I will model the beep to factor circadian rhythms
# Sex + physical_excercise_dur + acc_delta +
# temp_mean_z + temp_slope_z +
# (0+event_tot_s+activity_tot_s+social_tot_s+physical_tot_s | castor_record_id) +
# (0+ema_beep|castor_record_id)+
# (0+temp_slope_z|castor_record_id) + (0+temp_mean_z|castor_record_id)+
# (0+acc_delta|castor_record_id) + (0+physical_excercise_dur|castor_record_id),
# EMA_Data,
# family=Gamma(link="identity"),
# control=glmerControl(calc.derivs = FALSE) )
# Model Summary
model_output(glmer.sc_ph_auc_ema)
[4mMODEL INFO:[24m
[3mObservations:[23m 4565
[3mDependent Variable:[23m sc_phasic_auc_s
[3mType:[23m Mixed effects generalized linear regression
[3mError Distribution: [23mGamma
[3mLink function: [23midentity
[4mMODEL FIT:[24m
[3mAIC[23m = 9873.846, [3mBIC[23m = 10060.205
[3mPseudo-R² (fixed effects)[23m = NA
[3mPseudo-R² (total)[23m = NA
[4mFIXED EFFECTS:
[24m-------------------------------------------------------------------
Est. S.E. t val. p
--------------------------------- -------- ------- -------- -------
(Intercept) 1.061 0.158 6.700 0.000
event_tot_s 0.019 0.024 0.781 0.435
activity_tot_s 0.031 0.031 0.977 0.328
social_tot_s 0.014 0.022 0.644 0.520
physical_tot_s 0.011 0.011 0.998 0.318
ema_beep 0.010 0.009 1.085 0.278
SexMale 0.029 0.064 0.458 0.647
physical_excercise_dur 0.225 0.115 1.959 0.050
acc_delta 0.119 0.007 16.349 0.000
temp_mean_z 0.120 0.026 4.714 0.000
temp_slope_z 0.074 0.024 3.043 0.002
event_tot_s:activity_tot_s -0.004 0.005 -0.778 0.437
activity_tot_s:social_tot_s -0.004 0.004 -0.941 0.347
-------------------------------------------------------------------
[4m
RANDOM EFFECTS:
[24m---------------------------------------------------------
Group Parameter Std. Dev.
-------------------- ------------------------ -----------
castor_record_id event_tot_s 0.000
castor_record_id activity_tot_s 0.053
castor_record_id social_tot_s 0.050
castor_record_id physical_tot_s 0.052
castor_record_id.1 ema_beep 0.045
castor_record_id.2 temp_slope_z 0.135
castor_record_id.3 temp_mean_z 0.166
castor_record_id.4 acc_delta 0.052
castor_record_id.5 physical_excercise_dur 0.485
Residual 0.504
---------------------------------------------------------
[4m
Grouping variables:
[24m-------------------------------------
Group # groups ICC
------------------ ---------- -------
castor_record_id 82 0.000
-------------------------------------
Model has interaction terms. VIFs might be inflated. You may check multicollinearity among predictors of a model without interaction terms.
[34m# Check for Multicollinearity
[39m
[32mLow Correlation
[39m Parameter VIF Increased SE
physical_tot_s 1.13 1.06
ema_beep 1.05 1.02
Sex 1.01 1.00
physical_excercise_dur 1.00 1.00
acc_delta 1.01 1.00
temp_mean_z 1.02 1.01
temp_slope_z 1.00 1.00
[33mModerate Correlation
[39m Parameter VIF Increased SE
event_tot_s 5.21 2.28
social_tot_s 5.46 2.34
[31mHigh Correlation
[39m Parameter VIF Increased SE
activity_tot_s 14.59 3.82
event_tot_s:activity_tot_s 20.25 4.50
activity_tot_s:social_tot_s 13.38 3.66
Finally, we also take a look at subjective stress associations with subjective stress. Here we do not see any striking relationships.
# Model
# glmer.hr_mean_ema <- glmer(hr_mean_s ~
# event_tot_s*activity_tot_s + activity_tot_s*social_tot_s + physical_tot_s +
# ema_beep + # I will model the beep to factor circadian rhythms
# Sex + physical_excercise_dur + acc_delta +
# temp_mean_z + temp_slope_z +
# (0+event_tot_s+activity_tot_s+social_tot_s+physical_tot_s | castor_record_id) +
# (0+ema_beep|castor_record_id)+
# (0+temp_slope_z|castor_record_id) + (0+temp_mean_z|castor_record_id)+
# (0+acc_delta|castor_record_id) + (0+physical_excercise_dur|castor_record_id),
# EMA_Data,
# family=Gamma(link="identity"),
# control=glmerControl(calc.derivs = FALSE, optCtrl=list(maxfun=100000)))
# Model Summary
model_output(glmer.hr_mean_ema)
[4mMODEL INFO:[24m
[3mObservations:[23m 4419
[3mDependent Variable:[23m hr_mean_s
[3mType:[23m Mixed effects generalized linear regression
[3mError Distribution: [23mGamma
[3mLink function: [23midentity
[4mMODEL FIT:[24m
[3mAIC[23m = 8209.278, [3mBIC[23m = 8394.694
[3mPseudo-R² (fixed effects)[23m = NA
[3mPseudo-R² (total)[23m = NA
[4mFIXED EFFECTS:
[24m-------------------------------------------------------------------
Est. S.E. t val. p
--------------------------------- -------- ------- -------- -------
(Intercept) 3.205 0.120 26.733 0.000
event_tot_s 0.021 0.019 1.111 0.267
activity_tot_s -0.006 0.024 -0.271 0.786
social_tot_s -0.025 0.017 -1.456 0.145
physical_tot_s 0.023 0.010 2.257 0.024
ema_beep -0.042 0.007 -5.653 0.000
SexMale 0.020 0.064 0.313 0.754
physical_excercise_dur 0.359 0.080 4.482 0.000
acc_delta 0.144 0.005 29.562 0.000
temp_mean_z -0.089 0.022 -4.066 0.000
temp_slope_z 0.032 0.021 1.556 0.120
event_tot_s:activity_tot_s -0.003 0.004 -0.724 0.469
activity_tot_s:social_tot_s 0.004 0.003 1.210 0.226
-------------------------------------------------------------------
[4m
RANDOM EFFECTS:
[24m---------------------------------------------------------
Group Parameter Std. Dev.
-------------------- ------------------------ -----------
castor_record_id event_tot_s 0.052
castor_record_id activity_tot_s 0.040
castor_record_id social_tot_s 0.048
castor_record_id physical_tot_s 0.058
castor_record_id.1 ema_beep 0.039
castor_record_id.2 temp_slope_z 0.115
castor_record_id.3 temp_mean_z 0.147
castor_record_id.4 acc_delta 0.033
castor_record_id.5 physical_excercise_dur 0.282
Residual 0.172
---------------------------------------------------------
[4m
Grouping variables:
[24m-------------------------------------
Group # groups ICC
------------------ ---------- -------
castor_record_id 82 0.018
-------------------------------------
Model has interaction terms. VIFs might be inflated. You may check multicollinearity among predictors of a model without interaction terms.
[34m# Check for Multicollinearity
[39m
[32mLow Correlation
[39m Parameter VIF Increased SE
event_tot_s 4.51 2.12
social_tot_s 4.37 2.09
physical_tot_s 1.10 1.05
ema_beep 1.05 1.02
Sex 1.00 1.00
physical_excercise_dur 1.00 1.00
acc_delta 1.01 1.01
temp_mean_z 1.03 1.01
temp_slope_z 1.01 1.00
[31mHigh Correlation
[39m Parameter VIF Increased SE
activity_tot_s 13.76 3.71
event_tot_s:activity_tot_s 20.14 4.49
activity_tot_s:social_tot_s 10.72 3.27
# Model
# glmer.hr_min_ema <- glmer(hr_min_s ~
# event_tot_s*activity_tot_s + activity_tot_s*social_tot_s + physical_tot_s +
# ema_beep + # I will model the beep to factor circadian rhythms
# Sex + physical_excercise_dur + acc_delta +
# temp_mean_z + temp_slope_z +
# (0+event_tot_s+activity_tot_s+social_tot_s+physical_tot_s | castor_record_id) +
# (0+ema_beep|castor_record_id)+
# (0+temp_slope_z|castor_record_id) + (0+temp_mean_z|castor_record_id)+
# (0+acc_delta|castor_record_id) + (0+physical_excercise_dur|castor_record_id),
# EMA_Data,
# family=Gamma(link="log"),
# control=glmerControl(calc.derivs = FALSE) )
# Model Summary
model_output(glmer.hr_min_ema)
[4mMODEL INFO:[24m
[3mObservations:[23m 4419
[3mDependent Variable:[23m hr_min_s
[3mType:[23m Mixed effects generalized linear regression
[3mError Distribution: [23mGamma
[3mLink function: [23mlog
[4mMODEL FIT:[24m
[3mAIC[23m = 10278.955, [3mBIC[23m = 10464.372
[3mPseudo-R² (fixed effects)[23m = NA
[3mPseudo-R² (total)[23m = NA
[4mFIXED EFFECTS:
[24m-------------------------------------------------------------------
Est. S.E. t val. p
--------------------------------- -------- ------- -------- -------
(Intercept) 1.337 0.041 32.424 0.000
event_tot_s 0.003 0.007 0.479 0.632
activity_tot_s -0.009 0.008 -1.086 0.277
social_tot_s -0.016 0.006 -2.687 0.007
physical_tot_s 0.003 0.003 0.788 0.431
ema_beep -0.011 0.002 -4.569 0.000
SexMale -0.073 0.020 -3.576 0.000
physical_excercise_dur 0.109 0.030 3.660 0.000
acc_delta 0.022 0.001 19.155 0.000
temp_mean_z 0.016 0.007 2.256 0.024
temp_slope_z 0.009 0.006 1.422 0.155
event_tot_s:activity_tot_s -0.001 0.001 -0.787 0.431
activity_tot_s:social_tot_s 0.002 0.001 2.188 0.029
-------------------------------------------------------------------
[4m
RANDOM EFFECTS:
[24m---------------------------------------------------------
Group Parameter Std. Dev.
-------------------- ------------------------ -----------
castor_record_id event_tot_s 0.014
castor_record_id activity_tot_s 0.012
castor_record_id social_tot_s 0.013
castor_record_id physical_tot_s 0.021
castor_record_id.1 ema_beep 0.013
castor_record_id.2 temp_slope_z 0.035
castor_record_id.3 temp_mean_z 0.045
castor_record_id.4 acc_delta 0.007
castor_record_id.5 physical_excercise_dur 0.135
Residual 0.208
---------------------------------------------------------
[4m
Grouping variables:
[24m-------------------------------------
Group # groups ICC
------------------ ---------- -------
castor_record_id 82 0.003
-------------------------------------
Model has interaction terms. VIFs might be inflated. You may check multicollinearity among predictors of a model without interaction terms.
[34m# Check for Multicollinearity
[39m
[32mLow Correlation
[39m Parameter VIF Increased SE
event_tot_s 4.58 2.14
social_tot_s 4.91 2.22
physical_tot_s 1.12 1.06
ema_beep 1.05 1.02
Sex 1.00 1.00
physical_excercise_dur 1.00 1.00
acc_delta 1.02 1.01
temp_mean_z 1.03 1.01
temp_slope_z 1.01 1.00
[31mHigh Correlation
[39m Parameter VIF Increased SE
activity_tot_s 14.83 3.85
event_tot_s:activity_tot_s 20.81 4.56
activity_tot_s:social_tot_s 11.81 3.44
# Model
# glmer.hr_max_ema <- glmer(hr_max_s ~
# event_tot_s*activity_tot_s + activity_tot_s*social_tot_s + physical_tot_s +
# ema_beep + # I will model the beep to factor circadian rhythms
# Sex + physical_excercise_dur + acc_delta +
# temp_mean_z + temp_slope_z +
# (0+event_tot_s+activity_tot_s+social_tot_s+physical_tot_s | castor_record_id) +
# (0+ema_beep|castor_record_id)+
# (0+temp_slope_z|castor_record_id) + (0+temp_mean_z|castor_record_id)+
# (0+acc_delta|castor_record_id) + (0+physical_excercise_dur|castor_record_id),
# EMA_Data,
# family=Gamma(link="identity"),
# control=glmerControl(calc.derivs = FALSE) )
# Model Summary
model_output(glmer.hr_max_ema)
[4mMODEL INFO:[24m
[3mObservations:[23m 4419
[3mDependent Variable:[23m hr_max_s
[3mType:[23m Mixed effects generalized linear regression
[3mError Distribution: [23mGamma
[3mLink function: [23midentity
[4mMODEL FIT:[24m
[3mAIC[23m = 11315.568, [3mBIC[23m = 11500.984
[3mPseudo-R² (fixed effects)[23m = NA
[3mPseudo-R² (total)[23m = NA
[4mFIXED EFFECTS:
[24m-------------------------------------------------------------------
Est. S.E. t val. p
--------------------------------- -------- ------- -------- -------
(Intercept) 3.392 0.166 20.383 0.000
event_tot_s 0.024 0.027 0.879 0.379
activity_tot_s 0.025 0.033 0.748 0.455
social_tot_s -0.007 0.024 -0.276 0.783
physical_tot_s 0.011 0.013 0.846 0.397
ema_beep -0.056 0.010 -5.797 0.000
SexMale 0.159 0.068 2.337 0.019
physical_excercise_dur 0.337 0.134 2.515 0.012
acc_delta 0.174 0.006 28.427 0.000
temp_mean_z -0.241 0.029 -8.302 0.000
temp_slope_z 0.024 0.030 0.811 0.417
event_tot_s:activity_tot_s -0.004 0.005 -0.830 0.406
activity_tot_s:social_tot_s 0.001 0.004 0.129 0.898
-------------------------------------------------------------------
[4m
RANDOM EFFECTS:
[24m---------------------------------------------------------
Group Parameter Std. Dev.
-------------------- ------------------------ -----------
castor_record_id event_tot_s 0.074
castor_record_id activity_tot_s 0.050
castor_record_id social_tot_s 0.067
castor_record_id physical_tot_s 0.066
castor_record_id.1 ema_beep 0.042
castor_record_id.2 temp_slope_z 0.169
castor_record_id.3 temp_mean_z 0.188
castor_record_id.4 acc_delta 0.036
castor_record_id.5 physical_excercise_dur 0.657
Residual 0.219
---------------------------------------------------------
[4m
Grouping variables:
[24m-------------------------------------
Group # groups ICC
------------------ ---------- -------
castor_record_id 82 0.010
-------------------------------------
Model has interaction terms. VIFs might be inflated. You may check multicollinearity among predictors of a model without interaction terms.
[34m# Check for Multicollinearity
[39m
[32mLow Correlation
[39m Parameter VIF Increased SE
event_tot_s 4.88 2.21
social_tot_s 4.54 2.13
physical_tot_s 1.14 1.07
ema_beep 1.06 1.03
Sex 1.01 1.00
physical_excercise_dur 1.00 1.00
acc_delta 1.02 1.01
temp_mean_z 1.04 1.02
temp_slope_z 1.01 1.00
[31mHigh Correlation
[39m Parameter VIF Increased SE
activity_tot_s 16.25 4.03
event_tot_s:activity_tot_s 23.81 4.88
activity_tot_s:social_tot_s 12.09 3.48
tab_model(glmer.sc_ton_ema, glmer.sc_ph_num_ema, glmer.sc_ph_mag_ema, glmer.sc_ph_auc_ema,
glmer.hr_mean_ema, glmer.hr_min_ema, glmer.hr_max_ema,
dv.labels=c("SC Tonic", "SC Number", "SC Magnitude", "SC AUC", "Mean", "Min", "Max"),
title="Table X. Momentary Physiology and Subjective Stress",
transform=NULL, p.adjust="fdr",
show.stat=TRUE,
show.se=TRUE) %>%
return() %$%
knitr %>%
asis_output()
| SC Tonic | SC Number | SC Magnitude | SC AUC | Mean | Min | Max | |||||||||||||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | Statistic | p | Log-Mean | std. Error | CI | Statistic | p | Estimates | std. Error | CI | Statistic | p | Estimates | std. Error | CI | Statistic | p | Estimates | std. Error | CI | Statistic | p | Estimates | std. Error | CI | Statistic | p | Estimates | std. Error | CI | Statistic | p |
| (Intercept) | -0.00 | 0.05 | -0.11 – 0.10 | -0.09 | 0.993 | 2.01 | 0.11 | 1.78 – 2.23 | 17.69 | <0.001 | 1.34 | 0.16 | 1.01 – 1.66 | 8.12 | <0.001 | 1.06 | 0.16 | 0.75 – 1.37 | 6.70 | <0.001 | 3.21 | 0.12 | 2.97 – 3.44 | 26.73 | <0.001 | 1.34 | 0.04 | 1.26 – 1.42 | 32.42 | <0.001 | 3.39 | 0.17 | 3.07 – 3.72 | 20.38 | <0.001 |
| event_tot_s | 0.01 | 0.01 | -0.00 – 0.03 | 1.50 | 0.289 | 0.04 | 0.03 | -0.02 – 0.10 | 1.35 | 0.230 | 0.07 | 0.03 | 0.02 – 0.12 | 2.67 | 0.016 | 0.02 | 0.02 | -0.03 – 0.07 | 0.78 | 0.516 | 0.02 | 0.02 | -0.02 – 0.06 | 1.11 | 0.347 | 0.00 | 0.01 | -0.01 – 0.02 | 0.48 | 0.632 | 0.02 | 0.03 | -0.03 – 0.08 | 0.88 | 0.537 |
| activity_tot_s | 0.01 | 0.01 | -0.01 – 0.04 | 1.29 | 0.366 | 0.08 | 0.04 | 0.01 – 0.16 | 2.20 | 0.045 | 0.07 | 0.03 | 0.01 – 0.13 | 2.16 | 0.046 | 0.03 | 0.03 | -0.03 – 0.09 | 0.98 | 0.501 | -0.01 | 0.02 | -0.05 – 0.04 | -0.27 | 0.786 | -0.01 | 0.01 | -0.02 – 0.01 | -1.09 | 0.361 | 0.02 | 0.03 | -0.04 – 0.09 | 0.75 | 0.537 |
| social_tot_s | 0.00 | 0.01 | -0.02 – 0.02 | 0.01 | 0.993 | 0.03 | 0.03 | -0.03 – 0.09 | 0.90 | 0.396 | -0.00 | 0.02 | -0.05 – 0.04 | -0.18 | 0.931 | 0.01 | 0.02 | -0.03 – 0.06 | 0.64 | 0.563 | -0.03 | 0.02 | -0.06 – 0.01 | -1.46 | 0.236 | -0.02 | 0.01 | -0.03 – -0.00 | -2.69 | 0.016 | -0.01 | 0.02 | -0.05 – 0.04 | -0.28 | 0.848 |
| physical_tot_s | -0.00 | 0.00 | -0.01 – 0.01 | -0.53 | 0.869 | 0.04 | 0.03 | -0.02 – 0.09 | 1.16 | 0.293 | 0.03 | 0.01 | 0.00 – 0.05 | 2.14 | 0.046 | 0.01 | 0.01 | -0.01 – 0.03 | 1.00 | 0.501 | 0.02 | 0.01 | 0.00 – 0.04 | 2.26 | 0.052 | 0.00 | 0.00 | -0.00 – 0.01 | 0.79 | 0.467 | 0.01 | 0.01 | -0.01 – 0.04 | 0.85 | 0.537 |
| ema_beep | 0.00 | 0.00 | -0.00 – 0.01 | 0.43 | 0.869 | -0.05 | 0.02 | -0.10 – -0.00 | -2.14 | 0.047 | -0.00 | 0.01 | -0.02 – 0.02 | -0.07 | 0.944 | 0.01 | 0.01 | -0.01 – 0.03 | 1.09 | 0.501 | -0.04 | 0.01 | -0.06 – -0.03 | -5.65 | <0.001 | -0.01 | 0.00 | -0.02 – -0.01 | -4.57 | <0.001 | -0.06 | 0.01 | -0.08 – -0.04 | -5.80 | <0.001 |
| Sex [Male] | 0.00 | 0.02 | -0.04 – 0.04 | 0.06 | 0.993 | 0.84 | 0.11 | 0.62 – 1.07 | 7.36 | <0.001 | 0.04 | 0.07 | -0.10 – 0.18 | 0.58 | 0.728 | 0.03 | 0.06 | -0.10 – 0.15 | 0.46 | 0.647 | 0.02 | 0.06 | -0.11 – 0.15 | 0.31 | 0.786 | -0.07 | 0.02 | -0.11 – -0.03 | -3.58 | 0.001 | 0.16 | 0.07 | 0.03 – 0.29 | 2.34 | 0.042 |
| physical_excercise_dur | 0.08 | 0.05 | -0.01 – 0.18 | 1.69 | 0.261 | -0.20 | 0.41 | -1.00 – 0.60 | -0.48 | 0.630 | 0.36 | 0.13 | 0.11 – 0.62 | 2.79 | 0.014 | 0.22 | 0.11 | -0.00 – 0.45 | 1.96 | 0.130 | 0.36 | 0.08 | 0.20 – 0.52 | 4.48 | <0.001 | 0.11 | 0.03 | 0.05 – 0.17 | 3.66 | 0.001 | 0.34 | 0.13 | 0.07 – 0.60 | 2.52 | 0.031 |
| acc_delta | 0.02 | 0.00 | 0.01 – 0.02 | 9.52 | <0.001 | 0.23 | 0.01 | 0.20 – 0.26 | 16.04 | <0.001 | 0.11 | 0.01 | 0.10 – 0.12 | 18.21 | <0.001 | 0.12 | 0.01 | 0.10 – 0.13 | 16.35 | <0.001 | 0.14 | 0.00 | 0.13 – 0.15 | 29.56 | <0.001 | 0.02 | 0.00 | 0.02 – 0.02 | 19.15 | <0.001 | 0.17 | 0.01 | 0.16 – 0.19 | 28.43 | <0.001 |
| temp_mean_z | 0.04 | 0.01 | 0.02 – 0.06 | 4.49 | <0.001 | 0.43 | 0.09 | 0.25 – 0.61 | 4.71 | <0.001 | 0.12 | 0.02 | 0.07 – 0.17 | 4.81 | <0.001 | 0.12 | 0.03 | 0.07 – 0.17 | 4.71 | <0.001 | -0.09 | 0.02 | -0.13 – -0.05 | -4.07 | <0.001 | 0.02 | 0.01 | 0.00 – 0.03 | 2.26 | 0.045 | -0.24 | 0.03 | -0.30 – -0.18 | -8.30 | <0.001 |
| temp_slope_z | 0.03 | 0.01 | 0.01 – 0.05 | 2.84 | 0.020 | 0.22 | 0.07 | 0.09 – 0.36 | 3.24 | 0.003 | 0.07 | 0.02 | 0.03 – 0.11 | 3.23 | 0.004 | 0.07 | 0.02 | 0.03 – 0.12 | 3.04 | 0.008 | 0.03 | 0.02 | -0.01 – 0.07 | 1.56 | 0.222 | 0.01 | 0.01 | -0.00 – 0.02 | 1.42 | 0.224 | 0.02 | 0.03 | -0.03 – 0.08 | 0.81 | 0.537 |
|
event_tot_s * activity_tot_s |
-0.00 | 0.00 | -0.01 – 0.00 | -1.64 | 0.261 | -0.01 | 0.00 | -0.02 – -0.00 | -2.78 | 0.010 | -0.01 | 0.00 | -0.02 – -0.00 | -2.45 | 0.027 | -0.00 | 0.00 | -0.01 – 0.01 | -0.78 | 0.516 | -0.00 | 0.00 | -0.01 – 0.00 | -0.72 | 0.554 | -0.00 | 0.00 | -0.00 – 0.00 | -0.79 | 0.467 | -0.00 | 0.01 | -0.01 – 0.01 | -0.83 | 0.537 |
|
activity_tot_s * social_tot_s |
-0.00 | 0.00 | -0.00 – 0.00 | -0.45 | 0.869 | -0.01 | 0.00 | -0.02 – -0.01 | -3.92 | <0.001 | -0.00 | 0.00 | -0.01 – 0.01 | -0.32 | 0.883 | -0.00 | 0.00 | -0.01 – 0.00 | -0.94 | 0.501 | 0.00 | 0.00 | -0.00 – 0.01 | 1.21 | 0.327 | 0.00 | 0.00 | 0.00 – 0.00 | 2.19 | 0.047 | 0.00 | 0.00 | -0.01 – 0.01 | 0.13 | 0.898 |
| Random Effects | |||||||||||||||||||||||||||||||||||
| σ2 | 0.08 | 0.25 | 0.17 | 0.25 | 0.03 | 0.04 | 0.05 | ||||||||||||||||||||||||||||
| τ00 | 0.00 castor_record_id | 0.04 castor_record_id | 0.00 castor_record_id | 0.00 castor_record_id | 0.00 castor_record_id | 0.00 castor_record_id | 0.01 castor_record_id | ||||||||||||||||||||||||||||
| 0.00 castor_record_id.1 | 0.05 castor_record_id.1 | 0.00 castor_record_id.1 | 0.00 castor_record_id.1 | 0.00 castor_record_id.1 | 0.00 castor_record_id.1 | 0.00 castor_record_id.1 | |||||||||||||||||||||||||||||
| 0.00 castor_record_id.2 | 0.37 castor_record_id.2 | 0.01 castor_record_id.2 | 0.02 castor_record_id.2 | 0.01 castor_record_id.2 | 0.00 castor_record_id.2 | 0.03 castor_record_id.2 | |||||||||||||||||||||||||||||
| 0.00 castor_record_id.3 | 0.66 castor_record_id.3 | 0.02 castor_record_id.3 | 0.03 castor_record_id.3 | 0.02 castor_record_id.3 | 0.00 castor_record_id.3 | 0.04 castor_record_id.3 | |||||||||||||||||||||||||||||
| 0.00 castor_record_id.4 | 0.02 castor_record_id.4 | 0.00 castor_record_id.4 | 0.00 castor_record_id.4 | 0.00 castor_record_id.4 | 0.00 castor_record_id.4 | 0.00 castor_record_id.4 | |||||||||||||||||||||||||||||
| 0.08 castor_record_id.5 | 12.31 castor_record_id.5 | 0.39 castor_record_id.5 | 0.23 castor_record_id.5 | 0.08 castor_record_id.5 | 0.02 castor_record_id.5 | 0.43 castor_record_id.5 | |||||||||||||||||||||||||||||
| τ11 | 0.00 castor_record_id.activity_tot_s | 0.08 castor_record_id.activity_tot_s | 0.00 castor_record_id.activity_tot_s | 0.00 castor_record_id.activity_tot_s | 0.00 castor_record_id.activity_tot_s | 0.00 castor_record_id.activity_tot_s | 0.00 castor_record_id.activity_tot_s | ||||||||||||||||||||||||||||
| 0.00 castor_record_id.social_tot_s | 0.05 castor_record_id.social_tot_s | 0.00 castor_record_id.social_tot_s | 0.00 castor_record_id.social_tot_s | 0.00 castor_record_id.social_tot_s | 0.00 castor_record_id.social_tot_s | 0.00 castor_record_id.social_tot_s | |||||||||||||||||||||||||||||
| 0.00 castor_record_id.physical_tot_s | 0.07 castor_record_id.physical_tot_s | 0.00 castor_record_id.physical_tot_s | 0.00 castor_record_id.physical_tot_s | 0.00 castor_record_id.physical_tot_s | 0.00 castor_record_id.physical_tot_s | 0.00 castor_record_id.physical_tot_s | |||||||||||||||||||||||||||||
| ρ01 | -0.50 castor_record_id.activity_tot_s | -0.19 castor_record_id.activity_tot_s | -0.39 castor_record_id.activity_tot_s | -0.50 castor_record_id.activity_tot_s | -0.23 castor_record_id.activity_tot_s | -0.71 castor_record_id.activity_tot_s | |||||||||||||||||||||||||||||
| -0.58 castor_record_id.social_tot_s | -0.16 castor_record_id.social_tot_s | -0.39 castor_record_id.social_tot_s | -0.20 castor_record_id.social_tot_s | -0.14 castor_record_id.social_tot_s | -0.14 castor_record_id.social_tot_s | ||||||||||||||||||||||||||||||
| -0.54 castor_record_id.physical_tot_s | -0.23 castor_record_id.physical_tot_s | -0.26 castor_record_id.physical_tot_s | -0.38 castor_record_id.physical_tot_s | -0.36 castor_record_id.physical_tot_s | -0.62 castor_record_id.physical_tot_s | ||||||||||||||||||||||||||||||
| ICC | 0.07 | 0.89 | 0.33 | 0.72 | 0.15 | 0.59 | |||||||||||||||||||||||||||||
| N | 82 castor_record_id | 82 castor_record_id | 82 castor_record_id | 82 castor_record_id | 82 castor_record_id | 82 castor_record_id | 82 castor_record_id | ||||||||||||||||||||||||||||
| Observations | 4565 | 4565 | 4565 | 4565 | 4419 | 4419 | 4419 | ||||||||||||||||||||||||||||
| Marginal R2 / Conditional R2 | 0.102 / 0.169 | 0.409 / 0.935 | 0.533 / 0.688 | 0.577 / NA | 0.836 / 0.954 | 0.213 / 0.331 | 0.879 / 0.951 | ||||||||||||||||||||||||||||
Next I want to plot our outcomes as a function of the subjective stress. I will do this in separate graphs for event, activity, social, and physical stress. This gets a little complicated though. So I will make each graph separately, and then combine them outside of R to make sure they still look good.
# Shorrtlist the variables to those of interest
plot.vars <- plot_models(glmer.posmood_beep, glmer.sc_ph_num_ema, glmer.hr_mean_ema, glmer.hr_min_ema, glmer.hr_max_ema)
# Get and filter terms
full_term <- as.vector(plot.vars$data$term[1:length(plot.vars$data$term)])
event_term <- full_term[full_term != "event_tot_s"]
activ_term <-full_term[full_term != "activity_tot_s"]
soc_term <- full_term[full_term != "social_tot_s"]
phy_term <- full_term[full_term != "physical_tot_s"]
# Set the colours
col_list <-c('#7570B3','#7570B3','#7570B3',"#D95F02", "#D95F02", '#D95F02', '#D95F02', "#1B9E77", "#1B9E77")
# Event Plot
plot.event_substr <- plot_models(glmer.posmood_beep, glmer.negmood_beep,
glmer.sc_ton_ema, glmer.sc_ph_num_ema, glmer.sc_ph_mag_ema, glmer.sc_ph_auc_ema,
glmer.hr_mean_ema, glmer.hr_min_ema, glmer.hr_max_ema,
m.labels=c("Positive Affect", "Negative Affect", "SC Tonic", "SC Number", "SC Magnitued", "SC AUC", "HR Mean", "HR Minimum", "HR Maximum"),
axis.labels = c(""),
# Stastistical Stuff
rm.terms = event_term,
#show.values = T,
#value.size = 4,
#std.est=T,
show.p=T,
p.shape=T,
p.adjust = "fdr",
legend.pval.title = "Significance",
# Visual Stuff
colors=col_list,
dot.size=2,
line.size = 0.5,
spacing=0.7,
vline.color = "darkgrey",
legend.title = "") +
ptheme +
theme(axis.ticks.y=element_blank(),
axis.title = element_text(size=16),axis.text.x=element_text(size=11),
panel.grid.major.x = element_line(colour = "grey95"),
panel.grid.minor.x = element_line(colour = "grey90"))
plot.event_substr <- plot.event_substr+ coord_flip(ylim=c(0.75,1.25)) + coord_flip(ylim=c(-0.5,0.5)) + theme(panel.grid.major.y = element_blank(), panel.grid.minor.y = element_blank())
Coordinate system already present. Adding new coordinate system, which will replace the existing one.
Coordinate system already present. Adding new coordinate system, which will replace the existing one.
ggsave("figures/fig_StrResid_Event.tiff", units="in", width=4, height=4, dpi=300, compression = 'lzw', bg="transparent")
TIFFOpen: figures/fig_StrResid_Event.tiff: No such file or directory.
# Act Plot
plot.activity_substr <- plot_models(glmer.posmood_beep, glmer.negmood_beep,
glmer.sc_ton_ema, glmer.sc_ph_num_ema, glmer.sc_ph_mag_ema, glmer.sc_ph_auc_ema,
glmer.hr_mean_ema, glmer.hr_min_ema, glmer.hr_max_ema,
m.labels=c("Positive Affect", "Negative Affect", "SC Tonic", "SC Number", "SC Magnitued", "SC AUC", "HR Mean", "HR Minimum", "HR Maximum"),
axis.labels = c(""),
# Stastistical Stuff
rm.terms = activ_term,
#show.values = T,
#value.size = 4,
#std.est=T,
show.p=T,
p.shape=T,
p.adjust = "fdr",
legend.pval.title = "Significance",
# Visual Stuff
colors=col_list,
dot.size=2,
line.size = 0.5,
spacing=0.7,
vline.color = "darkgrey",
legend.title = "") +
ptheme +
theme(axis.ticks.y=element_blank(),
axis.title = element_text(size=16),axis.text.x=element_text(size=11),
panel.grid.major.x = element_line(colour = "grey95"),
panel.grid.minor.x = element_line(colour = "grey90"))
plot.activity_substr <- plot.activity_substr+ coord_flip(ylim=c(0.75,1.25)) + coord_flip(ylim=c(-0.5,0.5)) + theme(panel.grid.major.y = element_blank(), panel.grid.minor.y = element_blank())
ggsave("figures/fig_StrResid_Activity.tiff", units="in", width=4, height=4, dpi=300, compression = 'lzw', bg="transparent")
TIFFOpen: figures/fig_StrResid_Activity.tiff: No such file or directory.
# Social Stress
plot.social_substr <- plot_models(glmer.posmood_beep, glmer.negmood_beep,
glmer.sc_ton_ema, glmer.sc_ph_num_ema, glmer.sc_ph_mag_ema, glmer.sc_ph_auc_ema,
glmer.hr_mean_ema, glmer.hr_min_ema, glmer.hr_max_ema,
m.labels=c("Positive Affect", "Negative Affect", "SC Tonic", "SC Number", "SC Magnitued", "SC AUC", "HR Mean", "HR Minimum", "HR Maximum"),
axis.labels = c(""),
# Stastistical Stuff
rm.terms = soc_term,
#show.values = T,
#value.size = 4,
#std.est=T,
show.p=T,
p.shape=T,
p.adjust = "fdr",
legend.pval.title = "Significance",
# Visual Stuff
colors=col_list,
dot.size=2,
line.size = 0.5,
spacing=0.7,
vline.color = "darkgrey",
legend.title = "") +
ptheme +
theme(axis.ticks.y=element_blank(),
axis.title = element_text(size=16),axis.text.x=element_text(size=11),
panel.grid.major.x = element_line(colour = "grey95"),
panel.grid.minor.x = element_line(colour = "grey90"))
## Plot
plot.social_substr <- plot.social_substr + coord_flip(ylim=c(0.75,1.25)) + coord_flip(ylim=c(-0.5,0.5)) + theme(panel.grid.major.y = element_blank(), panel.grid.minor.y = element_blank())
ggsave("figures/fig_StrResid_Social.tiff", units="in", width=4, height=4, dpi=300, compression = 'lzw', bg="transparent")
TIFFOpen: figures/fig_StrResid_Social.tiff: No such file or directory.
# Physical Stress
plot.physical_substr <- plot_models(glmer.posmood_beep, glmer.negmood_beep,
glmer.sc_ton_ema, glmer.sc_ph_num_ema, glmer.sc_ph_mag_ema, glmer.sc_ph_auc_ema,
glmer.hr_mean_ema, glmer.hr_min_ema, glmer.hr_max_ema,
m.labels=c("Positive Affect", "Negative Affect", "SC Tonic", "SC Number", "SC Magnitued", "SC AUC", "HR Mean", "HR Minimum", "HR Maximum"),
axis.labels = c(""),
# Stastistical Stuff
rm.terms = phy_term,
#show.values = T,
#value.size = 4,
#std.est=T,
show.p=T,
p.shape=T,
p.adjust = "fdr",
legend.pval.title = "Significance",
# Visual Stuff
colors=col_list,
dot.size=2,
line.size =0.5,
spacing=0.7,
vline.color = "darkgrey",
legend.title = "") +
ptheme +
theme(axis.ticks.y=element_blank(),
axis.title = element_text(size=16),axis.text.x=element_text(size=11),
panel.grid.major.x = element_line(colour = "grey95"),
panel.grid.minor.x = element_line(colour = "grey90"))
# Plot
plot.physical_substr <- plot.physical_substr + coord_flip(ylim=c(-0.5,0.5), xlim=c(1,1.1)) + theme(panel.grid.major.y = element_blank(), panel.grid.minor.y = element_blank())
ggsave("figures/fig_StrResid_Physical.tiff", units="in", width=4, height=4, dpi=300, compression = 'lzw', bg="transparent")
TIFFOpen: figures/fig_StrResid_Physical.tiff: No such file or directory.
# All
ggarrange(plot.event_substr, plot.activity_substr, plot.social_substr, plot.physical_substr, common.legend = T, legend="right")
ggsave("figures/fig_StrResid_MAT.tiff", units="in", width=7, height=7, dpi=300, compression = 'lzw', bg="transparent")
TIFFOpen: figures/fig_StrResid_MAT.tiff: No such file or directory.
Arousal is not a valence specific mechanism. In the stress week, we see that there is a reduction in our physiological arousal measures, but also in positive affect. To this end, we check the relationship between physiology and the mood items. We expect that while arousal measures are indeed associated with stress, they are also associated with positive mood. To this end, we model these relationships with both positive and negative affect. We use a maximal fitting approach as before. We see a high correlation between positive and negative mood (presented below), and therefore also model interaction terms here. Due to model convergence issues, we only model the fixed intercepts here for the fixed effects of interest. Again, we check all model families before we pick the one with the best fit.
rcorr(as.matrix(EMA_Data[,c("mood_positive_s", "mood_negative_s")]))
mood_positive_s mood_negative_s
mood_positive_s 1.00 -0.57
mood_negative_s -0.57 1.00
n= 5735
P
mood_positive_s mood_negative_s
mood_positive_s 0
mood_negative_s 0
# Model
# glmer.sc_ton_mood <- glmer(sc_tonic_mean_s ~ mood_positive_s*mood_negative_s +
# ema_beep + # I will model the beep to factor circadian rhythms
# Sex + temp_slope_z + temp_mean_z +
# acc_delta + physical_excercise_dur +
# (0+mood_positive_s*mood_negative_s| castor_record_id) +# + (1+mood_negative_s|castor_record_id)+
# (0+ema_beep|castor_record_id)+
# (0+temp_slope_z|castor_record_id) + (0+temp_mean_z|castor_record_id)+
# (0+acc_delta|castor_record_id) + (0+physical_excercise_dur|castor_record_id),
# EMA_Data,
# family=Gamma(link="log"),
# control=glmerControl(calc.derivs = FALSE) )
# Summary
model_output(glmer.sc_ton_mood)
[4mMODEL INFO:[24m
[3mObservations:[23m 4565
[3mDependent Variable:[23m sc_tonic_mean_s
[3mType:[23m Mixed effects generalized linear regression
[3mError Distribution: [23mGamma
[3mLink function: [23mlog
[4mMODEL FIT:[24m
[3mAIC[23m = 773.006, [3mBIC[23m = 914.382
[3mPseudo-R² (fixed effects)[23m = NA
[3mPseudo-R² (total)[23m = NA
[4mFIXED EFFECTS:
[24m-----------------------------------------------------------------------
Est. S.E. t val. p
------------------------------------- -------- ------- -------- -------
(Intercept) 0.030 0.045 0.675 0.500
mood_positive_s 0.006 0.007 0.846 0.398
mood_negative_s -0.002 0.012 -0.145 0.885
ema_beep 0.001 0.003 0.252 0.801
SexMale -0.012 0.020 -0.637 0.524
temp_slope_z 0.029 0.010 2.934 0.003
temp_mean_z 0.042 0.009 4.532 0.000
acc_delta 0.019 0.002 9.752 0.000
physical_excercise_dur 0.082 0.050 1.631 0.103
mood_positive_s:mood_negative_s -0.001 0.002 -0.554 0.580
-----------------------------------------------------------------------
[4m
RANDOM EFFECTS:
[24m------------------------------------------------------------------
Group Parameter Std. Dev.
-------------------- --------------------------------- -----------
castor_record_id mood_positive_s 0.017
castor_record_id mood_negative_s 0.035
castor_record_id mood_positive_s:mood_negative_s 0.011
castor_record_id.1 ema_beep 0.016
castor_record_id.2 temp_slope_z 0.065
castor_record_id.3 temp_mean_z 0.062
castor_record_id.4 acc_delta 0.014
castor_record_id.5 physical_excercise_dur 0.299
Residual 0.287
------------------------------------------------------------------
[4m
Grouping variables:
[24m-------------------------------------
Group # groups ICC
------------------ ---------- -------
castor_record_id 82 0.002
-------------------------------------
[34m# Check for Multicollinearity
[39m
[32mLow Correlation
[39m Parameter VIF Increased SE
mood_positive_s 2.96 1.72
ema_beep 1.02 1.01
Sex 1.01 1.00
temp_slope_z 1.00 1.00
temp_mean_z 1.01 1.01
acc_delta 1.01 1.00
physical_excercise_dur 1.00 1.00
mood_positive_s:mood_negative_s 4.88 2.21
[33mModerate Correlation
[39m Parameter VIF Increased SE
mood_negative_s 6.44 2.54
# Model
# glmer.sc_phnum_mood <- glmer(sc_phasic_num_s ~ mood_positive_s*mood_negative_s +
# ema_beep + # I will model the beep to factor circadian rhythms
# Sex + temp_slope_z + temp_mean_z +
# acc_delta + physical_excercise_dur +
# (1+mood_positive_s*mood_negative_s| castor_record_id) +# + (1+mood_negative_s|castor_record_id)+
# (0+ema_beep|castor_record_id)+
# (0+temp_slope_z|castor_record_id) + (0+temp_mean_z|castor_record_id)+
# (0+acc_delta|castor_record_id) + (0+physical_excercise_dur|castor_record_id),
# EMA_Data,
# family=Gamma(link="identity"),
# control=glmerControl(calc.derivs = FALSE) )
# Summary
model_output(glmer.sc_phnum_mood)
[4mMODEL INFO:[24m
[3mObservations:[23m 4565
[3mDependent Variable:[23m log(sc_phasic_num_s)
[3mType:[23m Mixed effects linear regression
[4mMODEL FIT:[24m
[3mAIC[23m = 7580.862, [3mBIC[23m = 7747.943
[3mPseudo-R² (fixed effects)[23m = 0.202
[3mPseudo-R² (total)[23m = 0.331
[4mFIXED EFFECTS:
[24m---------------------------------------------------------------------------------
Est. S.E. t val. d.f. p
------------------------------------- -------- ------- -------- --------- -------
(Intercept) 0.232 0.086 2.690 205.786 0.008
mood_positive_s 0.029 0.011 2.555 343.443 0.011
mood_negative_s 0.037 0.020 1.842 54.171 0.071
ema_beep -0.017 0.005 -3.013 133.440 0.003
SexMale 0.073 0.057 1.273 67.946 0.207
temp_slope_z 0.036 0.009 4.016 5.496 0.008
temp_mean_z 0.074 0.015 4.874 63.065 0.000
acc_delta 0.056 0.003 17.256 85.710 0.000
physical_excercise_dur 0.084 0.061 1.378 64.266 0.173
mood_positive_s:mood_negative_s -0.006 0.004 -1.502 40.417 0.141
---------------------------------------------------------------------------------
[3m[36mp values calculated using Satterthwaite d.f.[39m[23m
[4m
RANDOM EFFECTS:
[24m------------------------------------------------------------------
Group Parameter Std. Dev.
-------------------- --------------------------------- -----------
castor_record_id (Intercept) 0.178
castor_record_id mood_positive_s 0.011
castor_record_id mood_negative_s 0.042
castor_record_id mood_positive_s:mood_negative_s 0.013
castor_record_id.1 ema_beep 0.022
castor_record_id.2 temp_slope_z 0.017
castor_record_id.3 temp_mean_z 0.091
castor_record_id.4 acc_delta 0.023
castor_record_id.5 physical_excercise_dur 0.162
Residual 0.523
------------------------------------------------------------------
[4m
Grouping variables:
[24m-------------------------------------
Group # groups ICC
------------------ ---------- -------
castor_record_id 82 0.093
-------------------------------------
[34m# Check for Multicollinearity
[39m
[32mLow Correlation
[39m Parameter VIF Increased SE
mood_positive_s 3.03 1.74
ema_beep 1.02 1.01
Sex 1.00 1.00
temp_slope_z 1.01 1.00
temp_mean_z 1.02 1.01
acc_delta 1.02 1.01
physical_excercise_dur 1.00 1.00
mood_positive_s:mood_negative_s 3.89 1.97
[33mModerate Correlation
[39m Parameter VIF Increased SE
mood_negative_s 5.76 2.40
# Model
# glmer.sc_ph_mag_mood <- glmer(sc_phasic_mag_s ~ mood_positive_s*mood_negative_s +
# ema_beep +
# Sex + temp_slope_z + temp_mean_z +
# acc_delta + physical_excercise_dur +
# (1+mood_positive_s*mood_negative_s| castor_record_id) +
# (0+ema_beep|castor_record_id)+
# (0+temp_slope_z|castor_record_id) + (0+temp_mean_z|castor_record_id)+
# (0+acc_delta|castor_record_id) + (0+physical_excercise_dur|castor_record_id),
# EMA_Data,
# family=Gamma(link="identity"),
# control=glmerControl(calc.derivs = FALSE, optCtrl=list(maxfun=100000) ))
# Summary
model_output(glmer.sc_ph_mag_mood)
[4mMODEL INFO:[24m
[3mObservations:[23m 4565
[3mDependent Variable:[23m sc_phasic_mag_s
[3mType:[23m Mixed effects generalized linear regression
[3mError Distribution: [23mGamma
[3mLink function: [23midentity
[4mMODEL FIT:[24m
[3mAIC[23m = 11181.350, [3mBIC[23m = 11348.431
[3mPseudo-R² (fixed effects)[23m = NA
[3mPseudo-R² (total)[23m = NA
[4mFIXED EFFECTS:
[24m-----------------------------------------------------------------------
Est. S.E. t val. p
------------------------------------- -------- ------- -------- -------
(Intercept) 1.571 0.160 9.802 0.000
mood_positive_s 0.035 0.022 1.640 0.101
mood_negative_s 0.029 0.035 0.831 0.406
ema_beep -0.006 0.009 -0.687 0.492
SexMale 0.020 0.077 0.259 0.796
temp_slope_z 0.070 0.021 3.355 0.001
temp_mean_z 0.126 0.024 5.217 0.000
acc_delta 0.104 0.006 17.239 0.000
physical_excercise_dur 0.348 0.127 2.731 0.006
mood_positive_s:mood_negative_s -0.006 0.007 -0.802 0.423
-----------------------------------------------------------------------
[4m
RANDOM EFFECTS:
[24m------------------------------------------------------------------
Group Parameter Std. Dev.
-------------------- --------------------------------- -----------
castor_record_id (Intercept) 0.679
castor_record_id mood_positive_s 0.084
castor_record_id mood_negative_s 0.116
castor_record_id mood_positive_s:mood_negative_s 0.035
castor_record_id.1 ema_beep 0.039
castor_record_id.2 temp_slope_z 0.086
castor_record_id.3 temp_mean_z 0.144
castor_record_id.4 acc_delta 0.039
castor_record_id.5 physical_excercise_dur 0.603
Residual 0.408
------------------------------------------------------------------
[4m
Grouping variables:
[24m-------------------------------------
Group # groups ICC
------------------ ---------- -------
castor_record_id 82 0.451
-------------------------------------
[34m# Check for Multicollinearity
[39m
[32mLow Correlation
[39m Parameter VIF Increased SE
mood_positive_s 3.37 1.84
ema_beep 1.02 1.01
Sex 1.01 1.00
temp_slope_z 1.01 1.00
temp_mean_z 1.02 1.01
acc_delta 1.01 1.01
physical_excercise_dur 1.00 1.00
mood_positive_s:mood_negative_s 3.37 1.83
[33mModerate Correlation
[39m Parameter VIF Increased SE
mood_negative_s 5.85 2.42
# Model
# glmer.sc_ph_auc_mood <- glmer(sc_phasic_auc_s ~ mood_positive_s*mood_negative_s +
# ema_beep +
# Sex + temp_slope_z + temp_mean_z +
# acc_delta + physical_excercise_dur +
# (1+mood_positive_s*mood_negative_s| castor_record_id) +
# (0+ema_beep|castor_record_id)+
# (0+temp_slope_z|castor_record_id) + (0+temp_mean_z|castor_record_id)+
# (0+acc_delta|castor_record_id) + (0+physical_excercise_dur|castor_record_id),
# EMA_Data,
# family=Gamma(link="log"),
# control=glmerControl(calc.derivs = FALSE) )
# Summary
model_output(glmer.sc_ph_auc_mood)
[4mMODEL INFO:[24m
[3mObservations:[23m 4565
[3mDependent Variable:[23m sc_phasic_auc_s
[3mType:[23m Mixed effects generalized linear regression
[3mError Distribution: [23mGamma
[3mLink function: [23mlog
[4mMODEL FIT:[24m
[3mAIC[23m = 9834.225, [3mBIC[23m = 10001.306
[3mPseudo-R² (fixed effects)[23m = NA
[3mPseudo-R² (total)[23m = NA
[4mFIXED EFFECTS:
[24m-----------------------------------------------------------------------
Est. S.E. t val. p
------------------------------------- -------- ------- -------- -------
(Intercept) 0.231 0.109 2.126 0.033
mood_positive_s 0.014 0.015 0.925 0.355
mood_negative_s 0.013 0.025 0.505 0.614
ema_beep -0.007 0.006 -1.183 0.237
SexMale 0.010 0.045 0.211 0.833
temp_slope_z 0.046 0.015 3.082 0.002
temp_mean_z 0.091 0.016 5.536 0.000
acc_delta 0.047 0.003 15.960 0.000
physical_excercise_dur 0.063 0.069 0.911 0.362
mood_positive_s:mood_negative_s -0.004 0.005 -0.779 0.436
-----------------------------------------------------------------------
[4m
RANDOM EFFECTS:
[24m------------------------------------------------------------------
Group Parameter Std. Dev.
-------------------- --------------------------------- -----------
castor_record_id (Intercept) 0.596
castor_record_id mood_positive_s 0.078
castor_record_id mood_negative_s 0.124
castor_record_id mood_positive_s:mood_negative_s 0.025
castor_record_id.1 ema_beep 0.029
castor_record_id.2 temp_slope_z 0.091
castor_record_id.3 temp_mean_z 0.111
castor_record_id.4 acc_delta 0.021
castor_record_id.5 physical_excercise_dur 0.312
Residual 0.492
------------------------------------------------------------------
[4m
Grouping variables:
[24m-------------------------------------
Group # groups ICC
------------------ ---------- -------
castor_record_id 82 0.496
-------------------------------------
[34m# Check for Multicollinearity
[39m
[32mLow Correlation
[39m Parameter VIF Increased SE
mood_positive_s 3.68 1.92
ema_beep 1.02 1.01
Sex 1.01 1.00
temp_slope_z 1.01 1.00
temp_mean_z 1.02 1.01
acc_delta 1.01 1.01
physical_excercise_dur 1.00 1.00
[33mModerate Correlation
[39m Parameter VIF Increased SE
mood_negative_s 8.71 2.95
mood_positive_s:mood_negative_s 5.05 2.25
# glmer.hr_mean_mood <- glmer(hr_mean_s ~ mood_positive_s*mood_negative_s +
# ema_beep + # I will model the beep to factor circadian rhythms
# Sex + temp_slope_z + temp_mean_z +
# acc_delta + physical_excercise_dur +
# (1+mood_positive_s*mood_negative_s| castor_record_id) + # (1+mood_positive_s| castor_record_id) + (1+mood_negative_s|castor_record_id)+
# (0+ema_beep|castor_record_id)+
# (0+temp_slope_z|castor_record_id) + (0+temp_mean_z|castor_record_id)+
# (0+acc_delta|castor_record_id) + (0+physical_excercise_dur|castor_record_id),
# EMA_Data,
# family=Gamma(link="log"),
# control=glmerControl(calc.derivs = FALSE) )
# Summary
model_output(glmer.hr_mean_mood )
[4mMODEL INFO:[24m
[3mObservations:[23m 4419
[3mDependent Variable:[23m hr_mean_s
[3mType:[23m Mixed effects generalized linear regression
[3mError Distribution: [23mGamma
[3mLink function: [23mlog
[4mMODEL FIT:[24m
[3mAIC[23m = 8547.250, [3mBIC[23m = 8713.485
[3mPseudo-R² (fixed effects)[23m = NA
[3mPseudo-R² (total)[23m = NA
[4mFIXED EFFECTS:
[24m-----------------------------------------------------------------------
Est. S.E. t val. p
------------------------------------- -------- ------- -------- -------
(Intercept) 1.104 0.039 28.317 0.000
mood_positive_s 0.015 0.005 2.830 0.005
mood_negative_s 0.012 0.009 1.346 0.178
ema_beep -0.016 0.002 -7.690 0.000
SexMale -0.055 0.020 -2.714 0.007
temp_slope_z 0.007 0.005 1.248 0.212
temp_mean_z -0.031 0.006 -4.811 0.000
acc_delta 0.033 0.001 26.656 0.000
physical_excercise_dur 0.092 0.025 3.657 0.000
mood_positive_s:mood_negative_s -0.002 0.001 -1.063 0.288
-----------------------------------------------------------------------
[4m
RANDOM EFFECTS:
[24m------------------------------------------------------------------
Group Parameter Std. Dev.
-------------------- --------------------------------- -----------
castor_record_id (Intercept) 0.209
castor_record_id mood_positive_s 0.026
castor_record_id mood_negative_s 0.040
castor_record_id mood_positive_s:mood_negative_s 0.007
castor_record_id.1 ema_beep 0.011
castor_record_id.2 temp_slope_z 0.028
castor_record_id.3 temp_mean_z 0.045
castor_record_id.4 acc_delta 0.009
castor_record_id.5 physical_excercise_dur 0.111
Residual 0.177
------------------------------------------------------------------
[4m
Grouping variables:
[24m-------------------------------------
Group # groups ICC
------------------ ---------- -------
castor_record_id 82 0.483
-------------------------------------
[34m# Check for Multicollinearity
[39m
[32mLow Correlation
[39m Parameter VIF Increased SE
mood_positive_s 3.36 1.83
ema_beep 1.02 1.01
Sex 1.01 1.00
temp_slope_z 1.01 1.00
temp_mean_z 1.02 1.01
acc_delta 1.01 1.01
physical_excercise_dur 1.00 1.00
[33mModerate Correlation
[39m Parameter VIF Increased SE
mood_negative_s 8.10 2.85
mood_positive_s:mood_negative_s 5.29 2.30
# glmer.hr_min_mood <- glmer(hr_min_s ~ mood_positive_s*mood_negative_s +
# ema_beep + # I will model the beep to factor circadian rhythms
# Sex + temp_slope_z + temp_mean_z +
# acc_delta + physical_excercise_dur +
# (1+mood_positive_s*mood_negative_s| castor_record_id) + # (1+mood_positive_s| castor_record_id) + (1+mood_negative_s|castor_record_id)+
# (0+ema_beep|castor_record_id)+
# (0+temp_slope_z|castor_record_id) + (0+temp_mean_z|castor_record_id)+
# (0+acc_delta|castor_record_id) + (0+physical_excercise_dur|castor_record_id),
# EMA_Data,
# family=Gamma(link="log"),
# control=glmerControl(calc.derivs = FALSE) )
# Summary
model_output(glmer.hr_min_mood )
[4mMODEL INFO:[24m
[3mObservations:[23m 4419
[3mDependent Variable:[23m hr_min_s
[3mType:[23m Mixed effects generalized linear regression
[3mError Distribution: [23mGamma
[3mLink function: [23mlog
[4mMODEL FIT:[24m
[3mAIC[23m = 10277.446, [3mBIC[23m = 10443.682
[3mPseudo-R² (fixed effects)[23m = NA
[3mPseudo-R² (total)[23m = NA
[4mFIXED EFFECTS:
[24m-----------------------------------------------------------------------
Est. S.E. t val. p
------------------------------------- -------- ------- -------- -------
(Intercept) 1.172 0.043 26.968 0.000
mood_positive_s 0.015 0.006 2.554 0.011
mood_negative_s 0.009 0.010 0.920 0.358
ema_beep -0.011 0.002 -4.680 0.000
SexMale -0.084 0.023 -3.591 0.000
temp_slope_z 0.010 0.006 1.503 0.133
temp_mean_z 0.020 0.007 2.782 0.005
acc_delta 0.022 0.001 18.734 0.000
physical_excercise_dur 0.105 0.030 3.540 0.000
mood_positive_s:mood_negative_s -0.001 0.002 -0.751 0.452
-----------------------------------------------------------------------
[4m
RANDOM EFFECTS:
[24m------------------------------------------------------------------
Group Parameter Std. Dev.
-------------------- --------------------------------- -----------
castor_record_id (Intercept) 0.214
castor_record_id mood_positive_s 0.027
castor_record_id mood_negative_s 0.043
castor_record_id mood_positive_s:mood_negative_s 0.007
castor_record_id.1 ema_beep 0.010
castor_record_id.2 temp_slope_z 0.035
castor_record_id.3 temp_mean_z 0.047
castor_record_id.4 acc_delta 0.007
castor_record_id.5 physical_excercise_dur 0.135
Residual 0.208
------------------------------------------------------------------
[4m
Grouping variables:
[24m-------------------------------------
Group # groups ICC
------------------ ---------- -------
castor_record_id 82 0.414
-------------------------------------
[34m# Check for Multicollinearity
[39m
[32mLow Correlation
[39m Parameter VIF Increased SE
mood_positive_s 3.12 1.77
ema_beep 1.03 1.01
Sex 1.01 1.00
temp_slope_z 1.01 1.00
temp_mean_z 1.02 1.01
acc_delta 1.02 1.01
physical_excercise_dur 1.00 1.00
mood_positive_s:mood_negative_s 4.84 2.20
[33mModerate Correlation
[39m Parameter VIF Increased SE
mood_negative_s 7.11 2.67
# glmer.hr_max_mood <- glmer(hr_max_s ~ mood_positive_s*mood_negative_s +
# ema_beep + # I will model the beep to factor circadian rhythms
# Sex + temp_slope_z + temp_mean_z +
# acc_delta + physical_excercise_dur +
# (1+mood_positive_s*mood_negative_s| castor_record_id) + # (1+mood_positive_s| castor_record_id) + (1+mood_negative_s|castor_record_id)+
# (0+ema_beep|castor_record_id)+
# (0+temp_slope_z|castor_record_id) + (0+temp_mean_z|castor_record_id)+
# (0+acc_delta|castor_record_id) + (0+physical_excercise_dur|castor_record_id),
# EMA_Data,
# family=Gamma(link="log"),
# control=glmerControl(calc.derivs = FALSE, #, optimizer="Nelder_Mead",
# optCtrl=list(maxfun=100000)) )
# Summary
model_output(glmer.hr_max_mood )
[4mMODEL INFO:[24m
[3mObservations:[23m 4419
[3mDependent Variable:[23m hr_max_s
[3mType:[23m Mixed effects generalized linear regression
[3mError Distribution: [23mGamma
[3mLink function: [23mlog
[4mMODEL FIT:[24m
[3mAIC[23m = 11629.004, [3mBIC[23m = 11795.240
[3mPseudo-R² (fixed effects)[23m = NA
[3mPseudo-R² (total)[23m = NA
[4mFIXED EFFECTS:
[24m-----------------------------------------------------------------------
Est. S.E. t val. p
------------------------------------- -------- ------- -------- -------
(Intercept) 1.197 0.043 27.672 0.000
mood_positive_s 0.015 0.006 2.557 0.011
mood_negative_s 0.014 0.009 1.621 0.105
ema_beep -0.017 0.003 -6.812 0.000
SexMale 0.019 0.017 1.093 0.275
temp_slope_z 0.003 0.007 0.434 0.664
temp_mean_z -0.065 0.007 -8.822 0.000
acc_delta 0.033 0.001 25.960 0.000
physical_excercise_dur 0.073 0.034 2.149 0.032
mood_positive_s:mood_negative_s -0.002 0.002 -1.405 0.160
-----------------------------------------------------------------------
[4m
RANDOM EFFECTS:
[24m------------------------------------------------------------------
Group Parameter Std. Dev.
-------------------- --------------------------------- -----------
castor_record_id (Intercept) 0.202
castor_record_id mood_positive_s 0.029
castor_record_id mood_negative_s 0.027
castor_record_id mood_positive_s:mood_negative_s 0.005
castor_record_id.1 ema_beep 0.011
castor_record_id.2 temp_slope_z 0.033
castor_record_id.3 temp_mean_z 0.049
castor_record_id.4 acc_delta 0.008
castor_record_id.5 physical_excercise_dur 0.169
Residual 0.226
------------------------------------------------------------------
[4m
Grouping variables:
[24m-------------------------------------
Group # groups ICC
------------------ ---------- -------
castor_record_id 82 0.330
-------------------------------------
[34m# Check for Multicollinearity
[39m
[32mLow Correlation
[39m Parameter VIF Increased SE
mood_positive_s 3.19 1.79
ema_beep 1.02 1.01
Sex 1.01 1.01
temp_slope_z 1.01 1.00
temp_mean_z 1.02 1.01
acc_delta 1.01 1.01
physical_excercise_dur 1.00 1.00
mood_positive_s:mood_negative_s 4.92 2.22
[33mModerate Correlation
[39m Parameter VIF Increased SE
mood_negative_s 7.64 2.76
tab_model(glmer.sc_ton_mood, glmer.sc_phnum_mood, glmer.sc_ph_mag_mood, glmer.sc_ph_auc_mood,
glmer.hr_mean_mood, glmer.hr_min_mood, glmer.hr_max_mood,
dv.labels=c("SC Tonic", "SC Number", "SC Magnitude", "SC AUC", "HR Mean", "HR Min", "HR Max"),
title="Table X. Physio vs. Mood",
transform=NULL, p.adjust="fdr",
show.stat=TRUE,
show.se=TRUE) %>%
return() %$%
knitr %>%
asis_output()
| SC Tonic | SC Number | SC Magnitude | SC AUC | HR Mean | HR Min | HR Max | |||||||||||||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | Statistic | p | Estimates | std. Error | CI | Statistic | p | Estimates | std. Error | CI | Statistic | p | Estimates | std. Error | CI | Statistic | p | Estimates | std. Error | CI | Statistic | p | Estimates | std. Error | CI | Statistic | p | Estimates | std. Error | CI | Statistic | p |
| (Intercept) | 0.03 | 0.04 | -0.06 – 0.12 | 0.67 | 0.724 | 0.23 | 0.09 | 0.06 – 0.40 | 2.69 | 0.014 | 1.57 | 0.16 | 1.26 – 1.88 | 9.80 | <0.001 | 0.23 | 0.11 | 0.02 – 0.44 | 2.13 | 0.084 | 1.10 | 0.04 | 1.03 – 1.18 | 28.32 | <0.001 | 1.17 | 0.04 | 1.09 – 1.26 | 26.97 | <0.001 | 1.20 | 0.04 | 1.11 – 1.28 | 27.67 | <0.001 |
| mood_positive_s | 0.01 | 0.01 | -0.01 – 0.02 | 0.85 | 0.724 | 0.03 | 0.01 | 0.01 – 0.05 | 2.55 | 0.018 | 0.04 | 0.02 | -0.01 – 0.08 | 1.64 | 0.168 | 0.01 | 0.01 | -0.02 – 0.04 | 0.93 | 0.518 | 0.01 | 0.01 | 0.00 – 0.02 | 2.83 | 0.008 | 0.01 | 0.01 | 0.00 – 0.03 | 2.55 | 0.015 | 0.02 | 0.01 | 0.00 – 0.03 | 2.56 | 0.021 |
| mood_negative_s | -0.00 | 0.01 | -0.02 – 0.02 | -0.14 | 0.885 | 0.04 | 0.02 | -0.00 – 0.08 | 1.84 | 0.094 | 0.03 | 0.04 | -0.04 – 0.10 | 0.83 | 0.528 | 0.01 | 0.03 | -0.04 – 0.06 | 0.50 | 0.682 | 0.01 | 0.01 | -0.01 – 0.03 | 1.35 | 0.223 | 0.01 | 0.01 | -0.01 – 0.03 | 0.92 | 0.397 | 0.01 | 0.01 | -0.00 – 0.03 | 1.62 | 0.150 |
| ema_beep | 0.00 | 0.00 | -0.01 – 0.01 | 0.25 | 0.885 | -0.02 | 0.01 | -0.03 – -0.01 | -3.01 | 0.006 | -0.01 | 0.01 | -0.02 – 0.01 | -0.69 | 0.547 | -0.01 | 0.01 | -0.02 – 0.00 | -1.18 | 0.473 | -0.02 | 0.00 | -0.02 – -0.01 | -7.69 | <0.001 | -0.01 | 0.00 | -0.02 – -0.01 | -4.68 | <0.001 | -0.02 | 0.00 | -0.02 – -0.01 | -6.81 | <0.001 |
| Sex [Male] | -0.01 | 0.02 | -0.05 – 0.03 | -0.64 | 0.724 | 0.07 | 0.06 | -0.04 – 0.18 | 1.27 | 0.203 | 0.02 | 0.08 | -0.13 – 0.17 | 0.26 | 0.796 | 0.01 | 0.05 | -0.08 – 0.10 | 0.21 | 0.833 | -0.05 | 0.02 | -0.09 – -0.02 | -2.71 | 0.010 | -0.08 | 0.02 | -0.13 – -0.04 | -3.59 | 0.001 | 0.02 | 0.02 | -0.02 – 0.05 | 1.09 | 0.305 |
| temp_slope_z | 0.03 | 0.01 | 0.01 – 0.05 | 2.93 | 0.011 | 0.04 | 0.01 | 0.02 – 0.05 | 4.02 | <0.001 | 0.07 | 0.02 | 0.03 – 0.11 | 3.35 | 0.002 | 0.05 | 0.01 | 0.02 – 0.08 | 3.08 | 0.007 | 0.01 | 0.01 | -0.00 – 0.02 | 1.25 | 0.236 | 0.01 | 0.01 | -0.00 – 0.02 | 1.50 | 0.166 | 0.00 | 0.01 | -0.01 – 0.02 | 0.43 | 0.664 |
| temp_mean_z | 0.04 | 0.01 | 0.02 – 0.06 | 4.53 | <0.001 | 0.07 | 0.02 | 0.04 – 0.10 | 4.87 | <0.001 | 0.13 | 0.02 | 0.08 – 0.17 | 5.22 | <0.001 | 0.09 | 0.02 | 0.06 – 0.12 | 5.54 | <0.001 | -0.03 | 0.01 | -0.04 – -0.02 | -4.81 | <0.001 | 0.02 | 0.01 | 0.01 – 0.03 | 2.78 | 0.009 | -0.07 | 0.01 | -0.08 – -0.05 | -8.82 | <0.001 |
| acc_delta | 0.02 | 0.00 | 0.02 – 0.02 | 9.75 | <0.001 | 0.06 | 0.00 | 0.05 – 0.06 | 17.26 | <0.001 | 0.10 | 0.01 | 0.09 – 0.12 | 17.24 | <0.001 | 0.05 | 0.00 | 0.04 – 0.05 | 15.96 | <0.001 | 0.03 | 0.00 | 0.03 – 0.04 | 26.66 | <0.001 | 0.02 | 0.00 | 0.02 – 0.02 | 18.73 | <0.001 | 0.03 | 0.00 | 0.03 – 0.04 | 25.96 | <0.001 |
| physical_excercise_dur | 0.08 | 0.05 | -0.02 – 0.18 | 1.63 | 0.257 | 0.08 | 0.06 | -0.04 – 0.20 | 1.38 | 0.187 | 0.35 | 0.13 | 0.10 – 0.60 | 2.73 | 0.013 | 0.06 | 0.07 | -0.07 – 0.20 | 0.91 | 0.518 | 0.09 | 0.03 | 0.04 – 0.14 | 3.66 | 0.001 | 0.11 | 0.03 | 0.05 – 0.16 | 3.54 | 0.001 | 0.07 | 0.03 | 0.01 – 0.14 | 2.15 | 0.053 |
|
mood_positive_s * mood_negative_s |
-0.00 | 0.00 | -0.01 – 0.00 | -0.55 | 0.724 | -0.01 | 0.00 | -0.01 – 0.00 | -1.50 | 0.166 | -0.01 | 0.01 | -0.02 – 0.01 | -0.80 | 0.528 | -0.00 | 0.00 | -0.01 – 0.01 | -0.78 | 0.545 | -0.00 | 0.00 | -0.00 – 0.00 | -1.06 | 0.288 | -0.00 | 0.00 | -0.00 – 0.00 | -0.75 | 0.452 | -0.00 | 0.00 | -0.01 – 0.00 | -1.41 | 0.200 |
| Random Effects | |||||||||||||||||||||||||||||||||||
| σ2 | 0.08 | 0.27 | 0.17 | 0.24 | 0.03 | 0.04 | 0.05 | ||||||||||||||||||||||||||||
| τ00 | 0.00 castor_record_id | 0.03 castor_record_id | 0.46 castor_record_id | 0.36 castor_record_id | 0.04 castor_record_id | 0.05 castor_record_id | 0.04 castor_record_id | ||||||||||||||||||||||||||||
| 0.00 castor_record_id.1 | 0.00 castor_record_id.1 | 0.00 castor_record_id.1 | 0.00 castor_record_id.1 | 0.00 castor_record_id.1 | 0.00 castor_record_id.1 | 0.00 castor_record_id.1 | |||||||||||||||||||||||||||||
| 0.00 castor_record_id.2 | 0.00 castor_record_id.2 | 0.01 castor_record_id.2 | 0.01 castor_record_id.2 | 0.00 castor_record_id.2 | 0.00 castor_record_id.2 | 0.00 castor_record_id.2 | |||||||||||||||||||||||||||||
| 0.00 castor_record_id.3 | 0.01 castor_record_id.3 | 0.02 castor_record_id.3 | 0.01 castor_record_id.3 | 0.00 castor_record_id.3 | 0.00 castor_record_id.3 | 0.00 castor_record_id.3 | |||||||||||||||||||||||||||||
| 0.00 castor_record_id.4 | 0.00 castor_record_id.4 | 0.00 castor_record_id.4 | 0.00 castor_record_id.4 | 0.00 castor_record_id.4 | 0.00 castor_record_id.4 | 0.00 castor_record_id.4 | |||||||||||||||||||||||||||||
| 0.09 castor_record_id.5 | 0.03 castor_record_id.5 | 0.36 castor_record_id.5 | 0.10 castor_record_id.5 | 0.01 castor_record_id.5 | 0.02 castor_record_id.5 | 0.03 castor_record_id.5 | |||||||||||||||||||||||||||||
| τ11 | 0.00 castor_record_id.mood_negative_s | 0.00 castor_record_id.mood_positive_s | 0.01 castor_record_id.mood_positive_s | 0.01 castor_record_id.mood_positive_s | 0.00 castor_record_id.mood_positive_s | 0.00 castor_record_id.mood_positive_s | 0.00 castor_record_id.mood_positive_s | ||||||||||||||||||||||||||||
| 0.00 castor_record_id.mood_positive_s:mood_negative_s | 0.00 castor_record_id.mood_negative_s | 0.01 castor_record_id.mood_negative_s | 0.02 castor_record_id.mood_negative_s | 0.00 castor_record_id.mood_negative_s | 0.00 castor_record_id.mood_negative_s | 0.00 castor_record_id.mood_negative_s | |||||||||||||||||||||||||||||
| 0.00 castor_record_id.mood_positive_s:mood_negative_s | 0.00 castor_record_id.mood_positive_s:mood_negative_s | 0.00 castor_record_id.mood_positive_s:mood_negative_s | 0.00 castor_record_id.mood_positive_s:mood_negative_s | 0.00 castor_record_id.mood_positive_s:mood_negative_s | 0.00 castor_record_id.mood_positive_s:mood_negative_s | ||||||||||||||||||||||||||||||
| ρ01 | 0.54 castor_record_id.mood_negative_s | 1.00 castor_record_id.mood_positive_s | -0.73 castor_record_id.mood_positive_s | -0.92 castor_record_id.mood_positive_s | -0.92 castor_record_id.mood_positive_s | -0.87 castor_record_id.mood_positive_s | -0.95 castor_record_id.mood_positive_s | ||||||||||||||||||||||||||||
| -0.85 castor_record_id.mood_positive_s:mood_negative_s | 0.59 castor_record_id.mood_negative_s | -0.56 castor_record_id.mood_negative_s | -0.77 castor_record_id.mood_negative_s | -0.79 castor_record_id.mood_negative_s | -0.67 castor_record_id.mood_negative_s | -0.87 castor_record_id.mood_negative_s | |||||||||||||||||||||||||||||
| -0.73 castor_record_id.mood_positive_s:mood_negative_s | -0.08 castor_record_id.mood_positive_s:mood_negative_s | 0.41 castor_record_id.mood_positive_s:mood_negative_s | 0.55 castor_record_id.mood_positive_s:mood_negative_s | 0.37 castor_record_id.mood_positive_s:mood_negative_s | 0.54 castor_record_id.mood_positive_s:mood_negative_s | ||||||||||||||||||||||||||||||
| ICC | 0.05 | 0.40 | 0.13 | 0.17 | 0.18 | 0.07 | |||||||||||||||||||||||||||||
| N | 82 castor_record_id | 82 castor_record_id | 82 castor_record_id | 82 castor_record_id | 82 castor_record_id | 82 castor_record_id | 82 castor_record_id | ||||||||||||||||||||||||||||
| Observations | 4565 | 4565 | 4565 | 4565 | 4419 | 4419 | 4419 | ||||||||||||||||||||||||||||
| Marginal R2 / Conditional R2 | 0.107 / 0.154 | 0.232 / NA | 0.499 / 0.699 | 0.178 / 0.284 | 0.450 / 0.543 | 0.219 / 0.359 | 0.390 / 0.432 | ||||||||||||||||||||||||||||
# Subset variables of interests
plot.physio <- plot_models(glmer.sc_ton_mood, glmer.sc_phnum_mood, glmer.hr_mean_mood)
rm_term <- as.vector(plot.physio$data$term[1:97])
rm_term <- rm_term[rm_term != "mood_positive_s" & rm_term!="mood_negative_s"]
# Set the colours
col_list <-c('#7570B3','#7570B3','#7570B3',"#D95F02", "#D95F02", '#D95F02', '#D95F02', "#1B9E77", "#1B9E77")
# Make the plot
plot.pos_mood_resid <- plot_models( glmer.sc_ton_mood, glmer.sc_phnum_mood, glmer.sc_ph_mag_mood, glmer.sc_ph_auc_mood,
glmer.hr_mean_mood, glmer.hr_min_mood, glmer.hr_max_mood,
m.labels=c("SC Tonic", "SC Number", "SC Magnitued", "SC AUC","HR Mean", "HR Minimum", "HR Maximum"),
axis.labels = c(" "),
# Stastistical Stuff
transform=NULL,
rm.terms = rm_term,
#show.values = T,
#value.size = 4,
#std.est=T,
show.p=T,
p.shape=T,
p.adjust = "fdr",
legend.pval.title = "Significance",
# Visual Stuff
colors=col_list,
dot.size=2,
line.size = 1,
spacing=0.7,
vline.color = "darkgrey",
legend.title = "") +
ptheme + scale_y_continuous(breaks=c( -0.2, -.1,0, .1,.2)) +
theme(axis.ticks.y=element_blank(),
axis.title = element_text(size=16),axis.text.x=element_text(size=11),
panel.grid.major.x = element_line(colour = "grey95"),
panel.grid.minor.x = element_line(colour = "grey90"))
Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.
# Plot and Save
plot.pos_mood_resid + coord_flip(ylim=c(-.25,.25)) + theme(panel.grid.major.y = element_line(colour = "grey95"), panel.grid.minor.y = element_line(colour = "grey90"))
Coordinate system already present. Adding new coordinate system, which will replace the existing one.
ggsave("figures/fig_MoodkResid_ALL.tiff", units="in", width=5, height=5, dpi=300, compression = 'lzw', bg="transparent")
TIFFOpen: figures/fig_MoodkResid_ALL.tiff: No such file or directory.
In the stress week, we see a decrease in both positive mood and physiology. However, we also see an increase in stress during this time too. This is confusing given the moment to moment relationship between subjective stress and physiology. One potential explanation is that the physiology is more strongly associate with positive arousal and opposed to negative arousal. So we need to actually confirm that what we see in the stress week (i.e. the physiology arousal decrease) is actually due to the positive mood changes. To this end, we perfrom a mediation analysis. We first need to filter out the nan’s here for some reason though. We cant apply this analysis to mixed level models, so we instead just check the simple models here.
library(mediation)
set.seed(123)
EMA_Sub <- EMA_Data %>% dplyr::select(castor_record_id, week_type, sc_phasic_mag_s, sc_phasic_num_s, mood_positive_s, event_tot_s, Sex, Program, # Model pop differences
First_Week, ema_beep_f, ema_day,
physical_excercise_dur, acc_delta,
temp_mean_z,temp_slope_z) %>% filter(complete.cases(.))
Two arousal measures were reduced in the stress week, but also associated with increased subjective stress. These measures were the number of skin conductance responses and the magnitude of skin conductance responses. So here we will check the mediating effect of positive mood.
First we check whether positive mood mediates the effects of week type of the number of responses. We see in the control week, there is a trend with positive affect mediating the number of skin conductace responses.
summary(results1)
Causal Mediation Analysis
Quasi-Bayesian Confidence Intervals
Mediator Groups: castor_record_id
Outcome Groups: castor_record_id
Output Based on Overall Averages Across Groups
Estimate 95% CI Lower 95% CI Upper p-value
ACME (control) -0.02325 -0.04792 0.00 0.057 .
ACME (treated) -0.00274 -0.02730 0.02 0.838
ADE (control) -0.29585 -0.39702 -0.20 <2e-16 ***
ADE (treated) -0.27533 -0.37472 -0.18 <2e-16 ***
Total Effect -0.29859 -0.39527 -0.20 <2e-16 ***
Prop. Mediated (control) 0.07693 -0.00190 0.18 0.057 .
Prop. Mediated (treated) 0.00920 -0.07380 0.10 0.838
ACME (average) -0.01300 -0.03232 0.01 0.175
ADE (average) -0.28559 -0.38601 -0.19 <2e-16 ***
Prop. Mediated (average) 0.04307 -0.01697 0.12 0.175
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Sample Size Used: 4565
Simulations: 5000
We next check the mediating effect of positive affect on skin coductance magnitde. Here we see that the magnitude of the responses is significanlty mediated by positive affect. It seems in both cases, these effects are being driven by the control week.
summary(results2)
Causal Mediation Analysis
Quasi-Bayesian Confidence Intervals
Mediator Groups: castor_record_id
Outcome Groups: castor_record_id
Output Based on Overall Averages Across Groups
Estimate 95% CI Lower 95% CI Upper p-value
ACME (control) -0.01727 -0.03291 0.00 0.032 *
ACME (treated) -0.00911 -0.02534 0.01 0.249
ADE (control) -0.17070 -0.23445 -0.11 <2e-16 ***
ADE (treated) -0.16253 -0.22615 -0.10 <2e-16 ***
Total Effect -0.17981 -0.24213 -0.12 <2e-16 ***
Prop. Mediated (control) 0.09571 0.01000 0.21 0.032 *
Prop. Mediated (treated) 0.04989 -0.03485 0.16 0.249
ACME (average) -0.01319 -0.02533 0.00 0.028 *
ADE (average) -0.16661 -0.22979 -0.10 <2e-16 ***
Prop. Mediated (average) 0.07280 0.00709 0.16 0.028 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Sample Size Used: 4565
Simulations: 5000
The results show that indeed, the results are mediated by positive mood in the stress week.
We were able to establish that physiology and mood both change in the weeks under chronic stress. That is good, because now we can see how well we can use this information in predictive models. To this end, we use the randomForest package plus some extra functions to test our prediction. The point of this analysis is to determine whether was can use the EPA/wristwatch measures on their own to classify which week participants are in, and to additionally check whether this would be as good as asking about mood. Finally, we will see whether the combination of both the mood and EPA measures is actually the best option.
Before we begin, we load the required library, and preset the seed to keep the analysis consistent. We additionally filter out the incomplete datapoints on only keep the full dataset
library(randomForest)
set.seed(123)
We first run the first analysis to establish the Out-Of-Bag (OOB) error rate without any cross validation. This is done to establish the models first, and fine tune the paramters (for example, the number of trees required for model stability and such). This is done with three models in mind:
Model 1. Mood Model 2. Physiology Model 3. Combination
Model 1 classifies week from the EMA mood outcomes, while Model 2 does the same using the EPA output. We then finally test these models using the combination thereof in Model 3. After we check these models, we can then go ahead and run the models with cross validation in the next sections.Based on the OOB, we see that the combined model does the best, while the physiology model does the worst.
This model classifies week type from the affect items
# Run random forest model for mood vars
forest.mood <- randomForest(Week_Type ~ mood_positive + mood_negative,
data=df.RandomForest,
ntree=5000,
importance = TRUE)
forest.mood
Call:
randomForest(formula = Week_Type ~ mood_positive + mood_negative, data = df.RandomForest, ntree = 5000, importance = TRUE)
Type of random forest: classification
Number of trees: 5000
No. of variables tried at each split: 1
OOB estimate of error rate: 43.45%
Confusion matrix:
Control Stress class.error
Control 1397 824 0.3710041
Stress 1128 1144 0.4964789
Now we can check what the physio data looks like.
# Run random forest full model
vars.forest.physio <- c( "sc_tonic_mean", "sc_phasic_mag", "sc_phasic_dur", "sc_phasic_auc", "sc_phasic_num","hr_mean", "hr_max", "hr_min", "hr_sd", "acc_delta", "temp_mean")
forest.full <- randomForest( as.formula(paste("Week_Type", "~", (paste(vars.forest.physio, collapse="+")))) ,
data=df.RandomForest,
ntree=5000,
importance = TRUE)
forest.full
Call:
randomForest(formula = as.formula(paste("Week_Type", "~", (paste(vars.forest.physio, collapse = "+")))), data = df.RandomForest, ntree = 5000, importance = TRUE)
Type of random forest: classification
Number of trees: 5000
No. of variables tried at each split: 3
OOB estimate of error rate: 46.16%
Confusion matrix:
Control Stress class.error
Control 1134 1087 0.4894192
Stress 987 1285 0.4344190
randomForest::importance(forest.full)
Control Stress MeanDecreaseAccuracy MeanDecreaseGini
sc_tonic_mean -7.3258933 20.849167 12.711316 237.6452
sc_phasic_mag -13.1342512 21.810438 8.633426 174.4823
sc_phasic_dur 22.6128487 -9.991554 18.481631 103.6931
sc_phasic_auc 0.5795613 13.941930 21.177497 229.3223
sc_phasic_num 12.7783302 3.930762 24.203246 124.0250
hr_mean 0.1237185 6.856160 7.719897 224.0363
hr_max -2.1026592 16.766449 16.354252 224.7909
hr_min -14.4347153 28.891901 13.847666 222.3630
hr_sd -2.0155031 9.912478 7.873558 220.4058
acc_delta 21.0259172 -7.252166 12.903214 239.5145
temp_mean 1.5120140 19.132939 14.799187 243.0713
Next up we check the entire data set and see if using both is good.
# Run the forest with most important features
forest.important <- randomForest(Week_Type ~ .,
data=df.RandomForest,
ntree=5000,
importance = TRUE)
forest.important
Call:
randomForest(formula = Week_Type ~ ., data = df.RandomForest, ntree = 5000, importance = TRUE)
Type of random forest: classification
Number of trees: 5000
No. of variables tried at each split: 3
OOB estimate of error rate: 40.71%
Confusion matrix:
Control Stress class.error
Control 1262 959 0.4317875
Stress 870 1402 0.3829225
In order to perform the validation, we try out a leave-one-beep-out (LOBO) analysis. This analysis will run a random forest model for each subject separately, leaving out one survey when running the model. We then test the model on the beep that was left out. We repeat this process for the subject until we’ve removed each survey/beep out once. This results in a series of predictions that give us an error level for each subject. In order to do this, we have developed a function that is provided in this github directory mentioned in the introducton. This will produce the classification errors for each subject.
We also need to test all our models against the true chance levels. One can assume in theory that the prediction error for a dichotomous variable is 50/50, but we dont like to do things the easy way. Instead, we will do 10000 iterations where we bootstrap resample the stress and control weeks. This will give us the real error rate with a sampling distribution. This was made into a function, but for some reason I couldnt get the function to run in parallel processing on our unix system. So instead I run the parallel processing of the bootstrap outside of a function. Its not pretty, but it works.
We also test the LOBO model using against the bootstrapped LOBO model. We also made a function to do this. The function estimates the p-value by comparing it to the permuation test distribution. The function also has the option to estimate the p-value from an assumed normal distribution, but we simply use the actual distribution for this. We additionally compute a combined p-value with this function using Stoufferś method and the poolr package.
This results in the three tabs for each model:
Some parts of the code are commented to make producing the notebook easier, we see that the mean LOBO error rate is 33.45 percent, and that the model performs relatively well on most subjects. Based on the estimated bootstrap error distribution of around 50%, we also show that this model performs above chance.
# Then with the full model
vars.forest.mood <- c("mood_positive_c", "mood_negative_c")
# forest.lobo.mood <- rf.lobo(Data=EMA_Data,
# SubNr="castor_record_id",
# xvars=vars.forest.mood,
# yvar="Week_Type",
# NoTree=5000)
as.data.frame(forest.lobo.mood$total$subject_errors)
psych::describe(forest.lobo.mood$total$subject_errors)
# BootstrapLOBO.
# { cl <- makeCluster(NoCores-1)
# registerDoSNOW(cl)
# iterations <- 10000
# print("Running, this will take a while...", quote=F)
# pb <- txtProgressBar(max = iterations, style = 3)
# progress <- function(n) setTxtProgressBar(pb, n)
# opts <- list(progress = progress)
# forest.boot.mood <- foreach(i=1:iterations, .options.snow = opts ) %dopar% {
# rf.BootstrapError(Data=EMA_Data,
# Shuffle="Week_Type",
# Iterations=1,
# SubjectId="castor_record_id",
# xvar=vars.forest.mood,
# yvar="Week_Type",
# Trees=500)
# }
# print("Done!", quote=F)
# stopCluster(cl)
# # Now merge them into one dataframe
# pb <- txtProgressBar(max =length(forest.boot.mood), style = 3)
# forest.boot.mood.merge <- as.data.frame(forest.boot.mood[[1]][2])
# for (i in 2:length(forest.boot.mood)) {
# setTxtProgressBar(pb, i)
# df.2merge <- as.data.frame(forest.boot.mood[[i]][2])
# forest.boot.mood.merge <- merge (df.2merge, forest.boot.mood.merge , by="DataFrame.id")
# }
# # After merging the results into a single data frame, lets also everage the error rate
# forest.boot.mood.merge$averageError <- rowMeans(forest.boot.mood.merge [, -which(names( forest.boot.mood.merge) %in% c("DataFrame.id"))], na.rm=T)
# saveRDS(forest.boot.mood.merge, file = "data/bootstrap_mood.rds")
# }
Now that we have run the bootstrap, we also save it just in case, and then we get the average error, and calcualte the true SD and SE from the sampling distribution.
# Print the Errors
boot.mood.descr <- psych::describe(forest.boot.mood.merge$averageError)
# Also get the true SD of bootstraps
boot.mood.descr$sd <- mean((psych::describe(forest.boot.mood.merge))$sd, na.rm=T)
NAs introduced by coercionno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf
boot.mood.descr$se <- mean((psych::describe(forest.boot.mood.merge))$se, na.rm=T)
NAs introduced by coercionno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf
boot.mood.descr
lobo.mood_null <- rf.null_test(forest.lobo.mood, forest.boot.mood.merge, model_type="LOBO", method="actual")
lobo.mood_null
$p_val
[1] 0.1597 0.0197 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
[22] 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
[43] 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
[64] 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 NA
$combined
combined p-values with: Stouffer's method
number of p-values combined: 82
test statistic: 33.193 ~ N(0,1)
adjustment: none
combined p-value: <2e-16
The second model I will test is whether we can predict week type from the physiology data alone. For this model, the LOBO error rate is 36.11, which is still better than chance level predictions, but worse than model 1.
# Then with the full model
# vars.forest.physio <- c("sc_tonic_mean", "sc_phasic_mag", "sc_phasic_dur", "sc_phasic_auc", "sc_phasic_num","hr_mean", "hr_max", "hr_min", "hr_sd", "acc_delta", "temp_mean")
# forest.lobo.physio <- rf.lobo(Data=EMA_Data,
# SubNr="castor_record_id",
# xvars=vars.forest.physio,
# yvar="Week_Type",
# NoTree=5000)
as.data.frame(forest.lobo.physio$total$subject_errors)
psych::describe(forest.lobo.physio$total$subject_errors)
We again do the bootstrap model, and then save it as an object for later use.
# # LOBO Bootstrapping
# { cl <- makeCluster(NoCores-1)
# registerDoSNOW(cl)
# iterations <- 10000
# print("Running, this will take a while...", quote=F)
# pb <- txtProgressBar(max = iterations, style = 3)
# progress <- function(n) setTxtProgressBar(pb, n)
# opts <- list(progress = progress)
# forest.boot.physio <- foreach(i=1:iterations, .options.snow = opts ) %dopar% {
# rf.BootstrapError(Data=EMA_Data,
# Shuffle="Week_Type",
# Iterations=1,
# SubjectId="castor_record_id",
# xvar=vars.forest.physio,
# yvar="Week_Type",
# Trees=500)
# }
# print("Done!", quote=F)
# stopCluster(cl)
#
# # Now merge them into one dataframe
# pb <- txtProgressBar(max =length(forest.boot.physio), style = 3)
# forest.boot.physio.merge <- as.data.frame(forest.boot.physio[[1]][2])
# for (i in 2:length(forest.boot.physio)) {
# setTxtProgressBar(pb, i)
# df.2merge <- as.data.frame(forest.boot.physio[[i]][2])
# forest.boot.physio.merge <- merge(df.2merge, forest.boot.physio.merge , by="DataFrame.id",
# suffixes=c((paste(i,".x")),(paste(i,".y"))))
# }
# # Merge it to mean
# forest.boot.physio.merge$averageError <- rowMeans(forest.boot.physio.merge [, -which(names( forest.boot.physio.merge) %in% c("DataFrame.id"))], na.rm=T)
# saveRDS( forest.boot.physio.merge, file = "data/bootstrap_physio.rds")
# }
We now calculate the bootstrap SE and SD
# Print the Errors
boot.physio.descr <- psych::describe(forest.boot.physio.merge$averageError)
# Also get the true SD of bootstraps
boot.physio.descr$sd <- mean((psych::describe(forest.boot.physio.merge))$sd, na.rm=T)
NAs introduced by coercionno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf
boot.physio.descr$se <- mean((psych::describe(forest.boot.physio.merge))$se, na.rm=T)
NAs introduced by coercionno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf
boot.physio.descr
So now that I have the averaged bootstrap error for the physio model, lets check if the LOBO is more than random with a test against the bootstrap estimates.
lobo.physio_null <- rf.null_test(forest.lobo.physio, forest.boot.physio.merge, model_type="LOBO", method="actual")
lobo.physio_null
$p_val
[1] 0.5632 0.0333 0.0001 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
[22] 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
[43] 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
[64] 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
$combined
combined p-values with: Stouffer's method
number of p-values combined: 80
test statistic: 32.62 ~ N(0,1)
adjustment: none
combined p-value: <2e-16
According to the t-test against the bootstrap error sample, the test is above chance level here too. I will again do a model comparison at the end of the model-testing sections.
Finally we check the full model using both physio and mood variables. This model perfroms better than the mood model. It also perfroms better than chance for each subject, with the pooled effect being significant.
# Then with the full model
# vars.forest.x <- c("mood_positive_c", "mood_negative_c", "sc_tonic_mean", "sc_phasic_mag", "sc_phasic_dur", "sc_phasic_auc", "sc_phasic_num","hr_mean", "hr_max", "hr_min", "hr_sd", "temp_mean")
# forest.lobo.full <- rf.lobo(Data=EMA_Data,
# SubNr="castor_record_id",
# xvars=vars.forest.x,
# yvar="Week_Type",
# NoTree=5000)
as.data.frame(forest.lobo.full$total$subject_errors)
psych::describe(forest.lobo.full$total$subject_errors)
Then I can check the bootstrap error for the model combined model as done before:
# {cl <- makeCluster(NoCores-1)
# registerDoSNOW(cl)
# iterations <- 10000
# print("Running, this will take a while...", quote=F)
# pb <- txtProgressBar(max = iterations, style = 3)
# progress <- function(n) setTxtProgressBar(pb, n)
# opts <- list(progress = progress)
# forest.boot.combi <- foreach(i=1:iterations, .options.snow = opts ) %dopar% {
# rf.BootstrapError(Data=EMA_Data,
# Shuffle="Week_Type",
# Iterations=1,
# SubjectId="castor_record_id",
# xvar=vars.forest.x,
# yvar="Week_Type",
# Trees=500)
# }
# print("Done!", quote=F)
# stopCluster(cl)
#
# # Now merge them into one dataframe
# pb <- txtProgressBar(max =length(forest.boot.combi), style = 3)
# forest.boot.combi.merge <- as.data.frame(forest.boot.physio[[1]][2])
# for (i in 2:length(forest.boot.combi)) {
# setTxtProgressBar(pb, i)
# df.2merge <- as.data.frame(forest.boot.combi[[i]][2])
# forest.boot.combi.merge <- merge(df.2merge, forest.boot.combi.merge , by="DataFrame.id",
# suffixes=c((paste(i,".x")),(paste(i,".y"))))
# }
# # Merge it to mean
# forest.boot.combi.merge$averageError <- rowMeans(forest.boot.combi.merge [, -which(names( forest.boot.combi.merge) %in% c("DataFrame.id"))], na.rm=T)
# saveRDS( forest.boot.combi.merge, file = "data/bootstrap_combi.rds")
# }
We then again calculate the SE and SD of the models
# Print the Errors
boot.combi.descr <- psych::describe(forest.boot.combi.merge$averageError)
# Also get the true SD of bootstraps
boot.combi.descr$sd <- mean((psych::describe(forest.boot.combi.merge))$sd, na.rm=T)
NAs introduced by coercionno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf
boot.combi.descr$se <- mean((psych::describe(forest.boot.combi.merge))$se, na.rm=T)
NAs introduced by coercionno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf
boot.combi.descr
And finally I can test the difference between the full model and the bootstrapped one. So we can see again, that our model only improves slightly with the addition of the physio variables (1%). That doesnt seem great. But we can also guess that within the control week there may be times when people are stressed, and vice versa for the stress week.We will do the null test again.
lobo.combi_null <- rf.null_test(forest.lobo.full, forest.boot.combi.merge, model_type="LOBO", method="actual")
lobo.combi_null
$p_val
[1] 0.1731 0.0037 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
[22] 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
[43] 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
[64] 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
$combined
combined p-values with: Stouffer's method
number of p-values combined: 80
test statistic: 32.837 ~ N(0,1)
adjustment: none
combined p-value: <2e-16
We now test the LOBO models against each other. We start with the highest error rates, and that is the mood and physiology models (models 1 and 2). We then test the mood against the combined LOBO model, Based on these results, we see that Model 3 > Model 1 > Model 2.
We filter the data and have to do a couple of manipulations to do a paired t-test.We see that the mood does significantly better
rf.result.df1 <- as.data.frame(forest.lobo.physio$total$subject_errors)
colnames(rf.result.df1)[1] <- "physio"
rf.result.df1$id <- rownames(rf.result.df1)
rf.result.df2 <- as.data.frame(forest.lobo.mood$total$subject_errors)
colnames(rf.result.df2)[1] <- "mood"
rf.result.df2$id <- rownames(rf.result.df2)
rf.result.merge <- merge(rf.result.df1, rf.result.df2, by="id")
tidy(t.test(x=rf.result.merge$physio, y=rf.result.merge$mood, paired=T))
Now I test the mood vs full model. It appears that the full model does a lot better than the one with mood only. So it is safe to say that Model 3 > Model 1 > Model 2
rf.result.df1 <- as.data.frame(forest.lobo.full$total$subject_errors)
colnames(rf.result.df1)[1] <- "full"
rf.result.df1$id <- rownames(rf.result.df1)
rf.result.df2 <- as.data.frame(forest.lobo.mood$total$subject_errors)
colnames(rf.result.df2)[1] <- "mood"
rf.result.df2$id <- rownames(rf.result.df2)
rf.result.merge <- merge(rf.result.df1, rf.result.df2, by="id")
tidy(t.test(x=rf.result.merge$full, y=rf.result.merge$mood, paired=T))
In order to determine whether the models can be generalized beyond the individial, we use a standard Leave-One-Subject-Out approach. In this analysis, we remove an entire subjects data set from the training set. We train our models on the remaining data, and then test the classification errors on the subject that was left out. This process is repeated until every subject has been removed once, thus providing error rates for each subject. For each model, we present the estimated LOSO error rates, and the test against the previously estimated null distribution.
Our first model looks at the mood ratings predicting which week participants are in. I estimate the LOSO error rates below. The model here doesnt do as well as the model with the lobo. Lets try to use the LOBO null distribution to see if its above chance levels next.
# First with mood items
# vars.forest.mood <- c("mood_positive", "mood_negative")
# forest.loso.mood <- rf.loso(Data=EMA_Data,
# SubNr="castor_record_id",
# xvars=vars.forest.mood,
# yvar="Week_Type",
# NoTree=5000)
as.data.frame(forest.loso.mood$error_rate)
psych::describe(forest.loso.mood$error_rate*100, na.rm=T)
Here I will try to check the LOSO model vs. the null distribution estimated from the bootstrap we did earler. The loso model did not perform better than the bootstrapped estimates.
forest.loso_mood_null <- rf.null_test(forest.loso.mood, forest.boot.mood.merge, model_type="LOSO", method="actual")
forest.loso_mood_null
$p_val
[1] 0.3563 0.5799 0.0000 0.7926 0.0000 0.0000 0.0000 0.0000 0.3076 0.8950 0.0012 0.0076 0.8176 0.5599 1.0000 0.0000 0.0177 0.0225 0.0833 0.0000 1.0000
[22] 0.0000 1.0000 0.0000 0.8348 0.7252 1.0000 0.0000 0.9897 0.0000 0.0004 0.0000 0.0000 0.1273 0.0001 0.9952 0.9466 0.0113 0.0000 0.9514 0.2565 0.0000
[43] 0.0000 0.0000 0.9906 0.0000 0.0000 0.0000 0.0000 1.0000 0.0000 0.0000 0.0000 1.0000 0.0022 0.0000 0.6286 0.9995 0.0000 1.0000 0.9948 0.0000 0.0020
[64] 1.0000 0.0001 1.0000 0.0016 0.2036 0.0000 0.0000 1.0000 0.9998 0.6477 1.0000 1.0000 0.0015 1.0000 0.0643 0.0000 0.0282 0.0006 0.0000 1.0000
$combined
combined p-values with: Stouffer's method
number of p-values combined: 83
test statistic: -Inf ~ N(0,1)
adjustment: none
combined p-value: 1
The second model is looking at the physio variables predicting week type. We see that it only achieves 47.83% success. We will test this now against the bootstrapped error rate. #### Model{-}
# Then with the full model
# vars.forest.physio <- c("sc_tonic_mean", "sc_phasic_mag", "sc_phasic_dur", "sc_phasic_auc", "sc_phasic_num","hr_mean", "hr_max", "hr_min", "hr_sd", "acc_delta", "temp_mean")
# forest.loso.physio <- rf.loso(Data=EMA_Data,
# SubNr="castor_record_id",
# xvars=vars.forest.physio,
# yvar="Week_Type",
# NoTree=5000)
as.data.frame(forest.loso.physio$error_rate)
There were 50 or more warnings (use warnings() to see the first 50)
psych::describe(forest.loso.physio$error_rate*100)
Next we check the physio variables against the null distribution. Here we see that it its also not signficantly better than chance.
forest.loso_physio_null <- rf.null_test(forest.loso.physio, forest.boot.physio.merge, model_type="LOSO", method="actual")
forest.loso_physio_null
$p_val
[1] 0.4703 0.8593 0.2199 0.4417 0.8771 0.0782 0.9687 0.4285 0.0932 NA 0.0021 0.1961 0.0539 0.7650 0.2076 0.0000 0.0000 1.0000 0.0000 0.9629 1.0000
[22] 0.0000 0.0000 0.0305 0.0000 0.2412 1.0000 0.3724 0.9777 0.0000 0.0090 1.0000 0.5665 0.0000 0.0073 0.9108 0.5631 0.9988 0.9466 0.0000 0.0006 0.1017
[43] 0.0000 0.9665 0.5722 0.3070 0.0000 1.0000 0.0000 0.0000 1.0000 0.0001 1.0000 0.9998 0.0000 0.1559 0.1032 1.0000 0.0233 1.0000 1.0000 0.0017 0.9999
[64] 0.5887 1.0000 1.0000 0.0244 0.0000 0.0000 0.9664 0.0005 0.0000 0.0166 1.0000 0.0106 1.0000 0.0000 0.9998 1.0000 0.9994
$combined
combined p-values with: Stouffer's method
number of p-values combined: 79
test statistic: -Inf ~ N(0,1)
adjustment: none
combined p-value: 1
FInally, the third model looks at everything (i.e. both physio and mood). This has an error rate of 43% which still sucks. Lets test this next against the bootstrap error
# Then with the full model
# vars.forest.x <- c("mood_positive_c", "mood_negative_c", "sc_tonic_mean", "sc_phasic_mag", "sc_phasic_dur", "sc_phasic_auc", "sc_phasic_num","hr_mean", "hr_max", "hr_min", "hr_sd", "acc_delta", "temp_mean")
# forest.loso.full <- rf.loso(Data=EMA_Data,
# SubNr="castor_record_id",
# xvars=vars.forest.x,
# yvar="Week_Type",
# NoTree=5000)
as.data.frame(forest.loso.full$error_rate)
psych::describe(forest.loso.full$error_rate*100)
Finally we test the model with both against the null distribution.This one is also not significantly better than chance.
forest.loso_full_null <- rf.null_test(forest.loso.full, forest.boot.combi.merge, model_type="LOSO", method="actual")
forest.loso_full_null
$p_val
[1] 0.0496 0.2528 0.0000 0.0485 0.0641 0.7740 0.0409 0.0000 0.0001 NA 0.0545 0.6610 0.6780 0.0003 1.0000 0.0000 0.0000 0.5450 0.5457 0.9957 0.0000
[22] 0.0000 0.0000 0.0000 0.0000 0.5653 0.8627 0.3768 0.0000 0.0000 0.0000 1.0000 0.0000 0.1334 0.0000 0.0000 0.0000 0.2496 0.0000 0.0000 0.0000 0.0000
[43] 0.0000 0.0062 0.0000 0.0000 0.0000 0.8257 0.0000 0.0000 0.3471 0.0000 0.0009 0.9658 0.0000 0.0002 0.0000 0.0000 0.0000 1.0000 0.0000 0.0000 0.0000
[64] 0.0002 0.0000 1.0000 0.0001 0.0001 0.0000 0.0000 0.0000 0.0629 0.0199 0.0000 1.0000 0.0155 0.0081 0.0000 0.6212 0.0078
$combined
combined p-values with: Stouffer's method
number of p-values combined: 79
test statistic: -Inf ~ N(0,1)
adjustment: none
combined p-value: 1
Based on all three tests for the models against the null we can conclude that the LOSO models do not in fact perform better than chance.
Based on the error rate, lets see if one model does a better job than the other:
rf.loso.mood <- as.data.frame(forest.loso.mood$error_rate)
There were 45 warnings (use warnings() to see them)
rf.loso.mood$id <- rownames(rf.loso.mood)
rf.loso.physio <- as.data.frame(forest.loso.physio$error_rate)
rf.loso.physio$id <- rownames(rf.loso.physio)
rf.loso.full <- as.data.frame(forest.loso.full$error_rate)
rf.loso.full$id <- rownames(rf.loso.full)
rf.result.merge <- merge(rf.loso.mood, rf.loso.physio, by="id")
rf.result.merge <- merge(rf.result.merge, rf.loso.full, by="id")
tidy(t.test(x=rf.result.merge$`forest.loso.mood$error_rate`, y=rf.result.merge$`forest.loso.physio$error_rate`, paired=T))
tidy(t.test(x=rf.result.merge$`forest.loso.mood$error_rate`, y=rf.result.merge$`forest.loso.full$error_rate`, paired=T))
tidy(t.test(x=rf.result.merge$`forest.loso.physio$error_rate`, y=rf.result.merge$`forest.loso.full$error_rate`, paired=T))
It looks like there is no significant difference between the error rates in the paired sample t-test in these three models. What could this mean? That regardless of what model we use at a population level, we still are not able to make decent predictions. I think that these models arent good at population level predictions, and this is pretty obnvious if we look at the within-subject variance from our linear models. Most of the variance comes from between subjects. Another approach is warranted.
We next waht to make a plot of my models with the standard error. We first merge them into one dataframe, then we can actually plot them.
# I will rename the bootstrap coloumns so I can add them to the full model dataframe for plotting
df_boot1 <- forest.boot.mood.merge[,c("DataFrame.id" ,"averageError")]
There were 18 warnings (use warnings() to see them)
df_boot1 <- df_boot1 %>% dplyr::rename(id=DataFrame.id, error=averageError)
df_boot1 <- tibble::add_column(df_boot1, model="boot.mood", .after="id")
df_boot2 <- forest.boot.physio.merge %>% dplyr::select(DataFrame.id, averageError) %>% dplyr::rename(id=DataFrame.id, error=averageError)
df_boot2 <- tibble::add_column(df_boot2, model="boot.physio", .after="id")
df_boot3 <- forest.boot.combi.merge %>% dplyr::select(DataFrame.id, averageError) %>% dplyr::rename(id=DataFrame.id, error=averageError)
df_boot3 <- tibble::add_column(df_boot3, model="boot.full", .after="id")
# Join the boot straps
df_boot <- full_join(df_boot1, df_boot2)
Joining, by = c("id", "model", "error")
df_boot <- full_join(df_boot, df_boot3)
Joining, by = c("id", "model", "error")
# Make the results and merge boot
df_loo <- gather(rf.combined.errors,model, error, loso.mood:lobo.full, factor_key=TRUE)
df_loo <- full_join(df_loo, df_boot)
Joining, by = c("id", "model", "error")
Column `model` joining factor and character vector, coercing into character vector
# Separate coloumns for plotting
df_loo <- separate(df_loo, model, into=c("Method", "Model"))
df_loo_sum <- data_summary(df_loo, "error", c("Method", "Model"))
# Set the SD and SE of bootstrap models to the correct ones of whole model
df_loo_sum$sd[df_loo_sum$Method=="boot" & df_loo_sum$Model=="mood"] <- boot.mood.descr$sd
df_loo_sum$se[df_loo_sum$Method=="boot" & df_loo_sum$Model=="mood"] <- boot.mood.descr$se
df_loo_sum$sd[df_loo_sum$Method=="boot" & df_loo_sum$Model=="physio"] <- boot.physio.descr$sd
df_loo_sum$se[df_loo_sum$Method=="boot" & df_loo_sum$Model=="physio"] <- boot.physio.descr$se
df_loo_sum$sd[df_loo_sum$Method=="boot" & df_loo_sum$Model=="full"] <- boot.combi.descr$sd
df_loo_sum$se[df_loo_sum$Method=="boot" & df_loo_sum$Model=="full"] <- boot.combi.descr$se
# Set the factor to make it pretty
df_loo_sum$Method <- factor(df_loo_sum$Method, levels =c("lobo", "loso", "boot"), labels = c("LOBO", "LOSO", "Bootstrap"))
df_loo_sum$Model<- factor(df_loo_sum$Model, levels =c("mood", "physio", "full"), labels = c("Model 1:\nMood", "Model 2:\n Physiology", "Model 3:\nCombination"))
figure.loo <- ggplot(df_loo_sum, aes(x=Model, y=error, colour=Method, fill=Method)) +
geom_bar(stat="summary", position = position_dodge2()) +
geom_errorbar(aes(ymin=error-se, ymax=error+se), width=.2, position=position_dodge(.9), colour="black") +
scale_y_continuous(breaks=c(0,10,20,30,40,50, 60, 70, 80), minor_breaks=c(5,15,25,35,45,55,65,75, 85) )+
scale_colour_brewer(palette="Set2") + scale_fill_brewer(palette="Set2") +
xlab("")+ ylab("Error Rate (%)\n") +
theme(text=element_text(size=16, family="Calibri"),
axis.line = element_line(size = 1, colour = "grey"),
panel.background = element_rect(fill="transparent"),
panel.grid.minor.y = element_line(colour="grey95"),
panel.grid.major.y = element_line(colour = "grey95"));
ggsave("figures/fig_RFModels.tiff", units="in", width=4, height=4, dpi=300, compression = 'lzw')
figure.loo
Now we want to see for each of the models, which one performs best. I first need to split the long dataframe I made for the figures into a wide format.
df_loo_wide <- pivot_wider(df_loo,id_cols = id, names_from=c(Method, Model), values_from=error)
Warning messages:
1: Unknown or uninitialised column: `completed`.
2: Unknown or uninitialised column: `completed`.
3: Unknown or uninitialised column: `completed`.
4: Unknown or uninitialised column: `completed`.
5: Unknown or uninitialised column: `completed`.
6: Unknown or uninitialised column: `completed`.
7: Unknown or uninitialised column: `completed`.
8: Unknown or uninitialised column: `completed`.
9: Unknown or uninitialised column: `completed`.
print(df_loo_wide)
write_csv(df_loo_wide, file = "data/LOO_DF.csv")
Error in write_csv(df_loo_wide, file = "data/LOO_DF.csv") :
unused argument (file = "data/LOO_DF.csv")
tidy(t.test(df_loo_wide$lobo_mood, df_loo_wide$boot_mood , paired = T))
tidy(t.test(df_loo_wide$loso_mood, df_loo_wide$boot_mood , paired = T))
tidy(t.test(df_loo_wide$lobo_mood, df_loo_wide$loso_mood , paired = T))
tidy(t.test(df_loo_wide$lobo_physio, df_loo_wide$boot_physio , paired = T))
tidy(t.test(df_loo_wide$loso_physio, df_loo_wide$boot_physio , paired = T))
tidy(t.test(df_loo_wide$lobo_physio, df_loo_wide$loso_physio , paired = T))
tidy(t.test(df_loo_wide$lobo_full, df_loo_wide$boot_full , paired = T))
tidy(t.test(df_loo_wide$loso_full, df_loo_wide$boot_full , paired = T))
tidy(t.test(df_loo_wide$lobo_full, df_loo_wide$loso_full , paired = T))
sessionInfo()
R version 3.6.1 (2019-07-05)
Platform: x86_64-conda_cos6-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)
Matrix products: default
BLAS/LAPACK: /home/cogaff/raytut/.conda/envs/r_env/lib/R/lib/libRblas.so
locale:
[1] en_US.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] poolr_0.8-2 extrafont_0.17 stargazer_5.2.2 factoextra_1.0.7 knitr_1.32 ggpubr_0.4.0 qgraph_1.6.5
[8] randomForest_4.6-14 doSNOW_1.0.18 snow_0.4-3 iterators_1.0.12 foreach_1.5.0 future_1.18.0 performance_0.6.1
[15] r2glmm_0.1.2 effectsize_0.3.1 sjPlot_2.8.7 jtools_2.1.0 olsrr_0.5.3 optimx_2020-4.2 ppcor_1.1
[22] MASS_7.3-51.6 corrr_0.4.2 RcmdrMisc_2.7-0 sandwich_2.5-1 car_3.0-8 carData_3.0-4 Hmisc_4.4-0
[29] Formula_1.2-3 survival_3.2-3 lattice_0.20-41 psycho_0.5.0 psych_1.9.12.31 boot_1.3-25 esmpack_0.1-14
[36] remotes_2.1.1 scales_1.1.1 zoo_1.8-8 data.table_1.12.8 janitor_2.0.1 hablar_0.3.0 broom_0.5.6
[43] plyr_1.8.6 forcats_0.5.0 stringr_1.4.0 dplyr_0.8.5 purrr_0.3.4 readr_1.3.1 tidyr_1.0.2
[50] tibble_3.0.1 ggplot2_3.3.3 tidyverse_1.3.0 readxl_1.3.1 lmerTest_3.1-2 lme4_1.1-23 Matrix_1.2-18
loaded via a namespace (and not attached):
[1] proto_1.0.0 tidyselect_1.0.0 htmlwidgets_1.5.1 grid_3.6.1 munsell_0.5.0 codetools_0.2-16
[7] statmod_1.4.34 withr_2.4.2 colorspace_2.0-0 rstudioapi_0.11 stats4_3.6.1 ggsignif_0.6.0
[13] Rttf2pt1_1.3.8 listenv_0.8.0 huge_1.3.4.1 emmeans_1.5.1 mnormt_1.5-6 farver_2.1.0
[19] coda_0.19-3 vctrs_0.3.1 generics_0.0.2 TH.data_1.0-10 clusterGeneration_1.3.4 xfun_0.22
[25] R6_2.5.0 arm_1.11-1 graphicalVAR_0.2.3 assertthat_0.2.1 multcomp_1.4-13 nnet_7.3-14
[31] texreg_1.37.5 gtable_0.3.0 globals_0.12.5 goftest_1.2-2 rlang_0.4.5 splines_3.6.1
[37] extrafontdb_1.0 rstatix_0.6.0 acepack_1.4.1 checkmate_2.0.0 yaml_2.2.1 reshape2_1.4.4
[43] abind_1.4-5 modelr_0.1.8 d3Network_0.5.2.1 backports_1.1.6 rsconnect_0.8.16 tools_3.6.1
[49] lavaan_0.6-6 ellipsis_0.3.1 RColorBrewer_1.1-2 gsubfn_0.7 Rcpp_1.0.4.6 base64enc_0.1-3
[55] rpart_4.1-15 pbapply_1.4-2 ggrepel_0.8.2 haven_2.3.1 cluster_2.1.0 fs_1.4.2
[61] magrittr_2.0.1 openxlsx_4.1.5 reprex_0.3.0 mvtnorm_1.1-1 whisker_0.4 sjmisc_2.8.6
[67] hms_0.5.3 evaluate_0.14 xtable_1.8-4 rio_0.5.16 sjstats_0.18.1 jpeg_0.1-8.1
[73] gridExtra_2.3 ggeffects_1.0.1 shape_1.4.5 compiler_3.6.1 crayon_1.4.1 minqa_1.2.4
[79] htmltools_0.5.0 corpcor_1.6.9 lubridate_1.7.9 DBI_1.1.0 sjlabelled_1.1.6 dbplyr_1.4.4
[85] cli_2.4.0 parallel_3.6.1 insight_0.13.2 igraph_1.2.5 BDgraph_2.62 pkgconfig_2.0.3
[91] numDeriv_2016.8-1.1 foreign_0.8-71 xml2_1.3.2 pbivnorm_0.6.0 estimability_1.3 rvest_0.3.5
[97] snakecase_0.11.0 digest_0.6.27 parameters_0.13.0 rmarkdown_2.7 cellranger_1.1.0 htmlTable_2.0.1
[103] nortest_1.0-4 curl_4.3 gtools_3.8.2 rjson_0.2.20 nloptr_1.2.2.2 lifecycle_0.2.0
[109] nlme_3.1-147 glasso_1.11 jsonlite_1.7.0 pillar_1.4.4 mlVAR_0.4.4 httr_1.4.1
[115] glue_1.4.0 bayestestR_0.9.0 zip_2.0.4 fdrtool_1.2.15 MplusAutomation_0.7-3 png_0.1-7
[121] pander_0.6.3 glmnet_4.0-2 class_7.3-17 stringi_1.4.6 blob_1.2.1 latticeExtra_0.6-29
[127] e1071_1.7-3
Warning messages:
1: Unknown or uninitialised column: `completed`.
2: Unknown or uninitialised column: `completed`.
3: Unknown or uninitialised column: `completed`.
4: Unknown or uninitialised column: `completed`.
5: Unknown or uninitialised column: `completed`.
6: Unknown or uninitialised column: `completed`.
7: Unknown or uninitialised column: `completed`.
8: Unknown or uninitialised column: `completed`.
9: Unknown or uninitialised column: `completed`.
Social Stress
Relates to stress derived from social interactions, or lack thereof.